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Abstract 

The  element-free  Galerkin  (EFG)  method  is  a  mesh-free  method  for  solv¬ 
ing  solid  mechanics  problems  with  an  approximation  based  only  on  nodes.  It 
is  developed  here  for  linear  elastic  fracture  problems.  Smoothing  of  mesh-free 
approximations  near  nonconvex  boundaries  is  done  by  three  methds:  (1)  the 
diffraction  method,  in  which  the  nodal  domain  of  influence  wrapped  a  short 
distance  around  a  boundary,  (2)  the  transparency  method,  which  is  described 
only  for  cracks,  yields  continuous  approximations  by  gradually  severing  the 
domains  of  influence  near  crack  tips,  and  (3)  the  “see-through”  method,  or 
continuous  line  criterion,  which  does  not  enforce  a  discontinuity  or  crack  if  the 
tip  is  within  the  domain  of  influence.  Two  methods  for  enriching  EFG  approxi¬ 
mations  for  linear  elastic  fracture  problems  are  described:  extrinsic  enrichment 
involves  adding  the  form  of  the  solution  to  the  trial  function;  for  intrinsic  en¬ 
richment,  the  EFG  basis  is  expanded  to  include  terms  from  the  near  tip  crack 
solution.  Several  problems  are  solved  to  illustrate  the  effectiveness  of  EFG 
crack  propagation  with  smoothing  and  enrichment. 


1  Introduction 

The  element-free  Galerkin  (EFG)  method  is  a  recently  developed  computational  tool 
which  has  been  used  extensively  for  computing  arbitrary  crack  paths.  EFG  uses  a 
Galerkin  scheme  for  approximating  the  solution  to  partial  differential  equations  with 
an  approximant  written  in  terms  of  a  set  of  nodes  and  the  surfaces  of  the  model.  This 
class  of  methods  is  often  called  meshless,  gridless  or  particle  methods  because  of  the 
absence  of  any  predefined  nodal  connectivity. 

Meshless  methods  have  been  in  existence  since  the  late  1970’s.  Lucy  (1977)  intro¬ 
duced  a  particle  method  called  smoothed  particle  hydrodynamics  (SPH)  for  modeling 
astrophysical  phenomena  and  Gingold  and  Monaghan  (1982)  and  Monaghan  (1982) 
used  this  method  for  problems  without  boundaries,  such  as  rotating  stars  and  dust 
clouds.  Libersky  and  Petschek  (1991)  extended  this  method  to  solve  solid  mechan¬ 
ics  problems.  Swegle,  Hicks,  and  Attaway  (1995)  noted  a  tensile  instability  existed 
in  SPH  and  Attaway,  Heinstein,  and  Swegle  (1994)  coupled  SPH  to  finite  elements 
through  a  contact  algorithm. 

A  separate  branch  of  meshless  methods  arose  from  the  work  of  Nayroles,  Touzot, 
and  Villon  (1992),  who  proposed  a  diffuse  element  method  (DEM)  using  a  basis 
function  and  a  weight  function  to  form  a  local  approximation  based  on  a  set  of  nodes. 
Belytschko,  Lu,  and  Gu  (1994)  recognized  that  this  approximation  was  actually  the 
moving  least  squares  (MLS)  approximation  of  Lancaster  and  Salkauskas  (1981)  and 
developed  a  similar  method  called  the  element-free  Galerkin  (EFG)  method  which  has 
been  used  extensively  for  fracture  and  crack  growth  (Belytschko,  Gu,  and  Lu,  1994; 
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Belytschko,  Lu,  Gu,  and  Tabbara,  1995;  Lu,  Belytschko,  and  Tabbara,  1995).  Liu, 
Jun,  and  Zhang  (1995)  also  proposed  a  meshless  method  called  the  reproducing  kernel 
particle  method  (RKPM)  with  an  approximation  based  on  a  convolution  integral.  The 
form  of  the  convolution  was  similar  to  SPH,  but  it  contains  a  correction  function  which 
corrects  for  consistency  near  boundaries  and  for  nonuniform  spacing.  Belytschko, 
Krongauz,  Organ,  Fleming,  and  Krysl  (1996)  have  shown  that  the  discrete  form  of 
the  convolution  integral  yields  approximants  which  are  identical  to  those  which  come 
from  MLS. 

Meshless  methods  based  on  partitions  of  unity  seem  to  provide  an  efficient  ve¬ 
hicle  for  performing  hp  adaptivity.  These  methods  include  /ip-clouds  (Duaxte  and 
Oden,  1996b)  and  the  paxtition  of  unity  finite  element  method  (PUFEM)  (Melenk 
and  Babuska,  1996).  The  hp- cloud  method  begins  with  a  partition  of  unity  based 
on  moving  least  squares  and  enhances  the  polynomial  order  of  the  approximation 
through  an  extrinsic  basis  which  can  be  added  locally  to  nodes. 

Other  meshless  methods  which  have  evolved  are  the  particle  in  cell  (PIC)  method 
(Sulsky,  Zhou,  and  Schreyer,  1995),  the  generalized  finite  difference  method  (Liszka 
and  Orkisz,  1980)  and  the  finite  point  method  (Onate,  Idelsohn,  Zienkiewicz,  Taylor, 
and  Sacco,  1996).  Belytschko,  Krongauz,  Organ,  Fleming,  and  Krysl  (1996)  provide 
a  comprehensive  review  of  the  state-of-the-art  in  meshless  methods  and  details  the 
relationship  between  several  of  the  key-methods. 

The  element-free  Galerkin  method  can  provide  an  excellent  complement  to  finite 
element  methods  in  situations  where  finite  elements  is  not  well  suited.  A  formula¬ 
tion  has  been  presented  in  Belytschko,  Organ,  and  Krongauz  (1995)  for  consistently 
coupling  EFG  and  FE  by  blending  the  approximations  in  a  transition  element.  This 
allows  t'he  speed  and  simplicity  of  finite  elements  to  be  exploited  where  possible  while 
allowing  EFG  to  be  used  in  regions  where  a  meshless  method  is  appropriate. 

One  class  of  problems  which  is  inherently  difficult  with  finite  element  methods 
is  crack  propagation  for  arbitrary  crack  paths.  The  predefined  connectivity  of  the 
element  structure  in  finite  element  methods  requires  special  treatment  if  a  crack  is 
to  be  extended.  Crack  propagation  has  been  modeled  with  finite  elements  in  several 
ways.  Early  crack  growth  with  finite  elements  was  done  by  nodal  release  methods  in 
which  cracks  were  grown  along  finite  element  boundaries  (Malluck  and  King,  1980). 
As  the  crack  grew,  the  elements  were  unzipped  to  create  new  surfaces.  This  method 
suffers  from  a  great  deal  of  mesh  dependence  for  arbitrary  crack  propagation.  In 
addition,  the  crack  increment  is  set  by  the  element  size  because  the  crack  grows  from 
node  to  node. 

The  most  obvious  way  to  handle  crack  propagation  is  by  remeshing  the  geometry. 
Swenson  and  Ingraffea  (1988)  presented  a  local  remeshing  technique  in  which  elements 
ahead  of  the  crack  tip  in  the  propagation  direction  were  removed  and  the  crack  was 
extended.  This  area  was  triangulated  to  create  a  new  local  crack  tip  mesh.  This 
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method  has  the  advantage  that  a  vast  amount  of  existing  finite  element  technology 
can  be  used  for  support.  There  are  some  drawbacks  to  this  method  such  as  difficulties 
if  the  crack  step  size  is  too  small,  then  remeshing  is  unwarranted.  Complex  geometries 
and  interacting  crack  tips  are  also  potential  difficulties  as  is  rerunning  a  step  with  a 
smaller  crack  increment. 

Al-Ostaz  and  Jasiuk  (1995)  modeled  fracture  and  crack  growth  with  finite  elements 
by  deleting  elements  which  met  a  yield  criterion.  This  approach,  which  is  not  based 
on  fracture  mechanics,  requires  a  very  fine  mesh  to  get  an  acceptable  representation  of 
a  crack.  Other  techniques  for  modeling  crack  growth  include  spring  network  models 
in  which  the  material  is  represented  by  a  network  of  springs  and  crack  propagation 
is  simulated  by  breaking  springs  (Schlangen  and  Mier,  1992).  Boundary  element 
methods  have  also  been  used  to  model,  crack  propagation  (Gallego  and  Dominguez, 
1992).  This  method  is  attractive  due  to  the  absence  of  a  domain  mesh,  making  crack 
extension  relatively  simple.  However,  the  requirements  of  a  Green’s  function  has 
limited  the  scope  of  problems  this  method  is  able  to  solve. 

This  report  focuses  on  the  element-free  Galerkin  methods  for  computational  frac¬ 
ture  mechanics.  In  Section  1,  the  moving  least  squares  methodology  is  reviewed  and 
used  to  construct  discrete  EFG  approximations.  Issues  in  modeling  cracks  in  a  mesh¬ 
less  method  are  discussed.  The  elastostatic  boundary  value  problem  is  presented 
along  with  its  associated  weak  form.  -Approximation  of  the  solution  with  EFG  is 
presented  and  the  topics  of  nodal  domains  of  influence  and  integration  of  the  weak 
form  are  discussed.  It  should  be  noted  that  the  definition  of  a  meshless  method  is  one 
in  which  no  predefined  element  connectivity  exists  for  determining  the  approximant. 
Some  confusion  invariably  arises  when  the  meshless  approximant  is  used  in  a  Galerkin 
method.  Integration  of  the  weak  form  is  performed  by  Gauss  quadrature  which  re¬ 
quires  some  sort  of  integration  cells.  Although  this  detracts  from  the  “meshless” 
nature  of  the  method,  the  background  cell  structure  by  no  means  destroys  it. 

Smoothing  of  EFG  approximation  near  nonconvex  boundaries  is  presented  in  Sec¬ 
tion  4.  Without  smoothing,  EFG  approximations  near  nonconvex  boundaries  such 
as  crack  tips  will  contain  discontinuities  which  extend  into  the  domain.  These  dis¬ 
continuities  arise  due  to  nodal  domains  of  influence  which  are  truncated  whenever  a 
ray  from  a  node  to  a  sampling  point  grazes  the  boundary.  The  diffraction  method 
smooths  EFG  approximations  by  wrapping  the  nodal  support  a  short  distance  around 
the  point  at  which  the  discontinuity  would  begin.  The  transparency  method,  which 
is  written  only  for  cracks,  yields  smooth  approximations  by  gradually  enforcing  the 
crack.  The  “see-through”  method,  or  continuous  line  criterion,  does  not  enforce  non¬ 
convex  boundaries  for  situations  in  which  a  continuous  line  can  be  drawn  between 
the  node  and  a  sampling  point  without  leaving  the  domain  of  influence. 

Section  5  presents  enrichment  techniques  for  the  EFG  method.  These  methods 
hinge  on  knowledge  of  the  solution  form  and  are  developed  for  linear  elastic  cracks. 
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The  enrichment  method  can  be  categorized  in  two  ways:  extrinsic  enrichment  en¬ 
hances  the  approximation  by  adding  functions  to  the  approximation;  in  intrinsic 
enrichment,  the  intrinsic  basis  is  expanded  to  include  enrichment  terms.  For  linear 
elastic  cracks,  enrichment  can  consist  of  adding  the  y/r  term  or  adding  the  entire 
near-tip  asymptotic  solution. 

Several  example  problems  axe  presented  in  Section  6  to  show  the  effects  of  the 
smoothing  and  enrichment  as  well  as  illustrate  the  effectiveness  of  EFG  for  solving 
problems  of  arbitrary  crack  growth. 

2  Meshless  Approximations  by  MLS 

A  meshless  approximation  for  a  discrete  system  is  one  which  is  written  entirely  in 
terms  of  the  parameter  values  at  nodes  -  no  predefined  connectivities  axe  used  as 
with  finite  elements  or  finite  differences.  A  smooth,  monotonically  decreasing  weight 
function  is  defined  at  each  node  such  that  the  whole  domain  is  covered.  Common 
shapes  of  weight  functions  in  two  dimensions  are  circles  and  rectangles  (see  Figs.  1). 
The  counterparts  of  these  in  three  dimensions  axe  spheres  and  bricks. 

Two  commonly  used  weight  functions  axe  the  Gaussian  and  the  quaxtic  spline 
given  in  Eqs.  (1).  For  the  circular  weights  shown  in  Fig.  la,  the  weight  functions  are 


Gaussian: 


quartic  spline: 


where  dj  =  ||x  —  X/||  is  the  distance  from  a  sampling  point  x  to  a  node  Xj,  and  dmi 
is  the  domain  of  influence  of  a  node,  i.e.,  the  area  for  which  the  weight  function  is 
nonzero.  The  variable  c  in  the  Gaussian  weight  is  used  to  control  the  dilation  of 
the  weight  function.  It  is  useful  to  define  a  characteristic  nodal  spacing,  cj,  which 
is  a  distance  such  that  a  node  will  yield  a  minimum  set  of  neighbors  sufficient  for 
regularity  of  the  equations  used  to  determine  the  approximant.  The  weight  function 
parameters  are  defined  in  terms  of  cj 


w(di)  = 


u>(d/)  = 


exp (— (d7/c)2)  -  exp (— (dm//c)2) 
1  -  exp(— (dm//c)2) 


dj  ^  dmj 
dj  dm  i 


(la) 


\dmI )  \dml  J  \dmlj 


dj  >  dmj 


(ib) 


dmi  —  dmaxCj 

c  =  acj 


(2) 

(3) 
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(a)  Circular. 


Figure  1:  A  computational  model  for  a  meshless  method  showing  the  bound 
ary,  nodes  and  supports. 
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where  dmax  and  a  are  constants.  For  the  Gaussian  weight  function  in  Eq.  (la),  the 
parameter  a  is  a  dilation  parameter  which  controls  the  shape  of  the  weight  function. 
If  a  is  kept  constant  while  dmax  is  increased,  the  shape  of  the  weight  function  will  not 
change  and  the  effective  domain  of  influence  will  be  smaller  than  the  actual  domain  of 
influence.  It  is  necessary  to  have  the  ratio  dmax/oi  >  4.0  to  avoid  poorly-formed  shape 
functions.  In  addition,  a  >  0.5  is  needed  for  smooth  shape  functions  and  derivatives. 
In  this  report,  the  Gaussian  weight  with  dmax  =  2.5,  a  =  0.625(=  dmax/4)is  used 
unless  otherwise  stated.  The  characteristic  nodal  spacing,  C/,  is  chosen  as  the  distance 
to  the  second  nearest  neighbor  for  regularly  spaced  nodes  and  the  distance  to  the  third 
nearest  neighbor  for  irregularly  spaced  nodes. 

The  approximation  of  a  function  u(x)  at  any  point  x  in  the  domain  Ll  is  written 

«*(x)  =  pT(x)a(x)  (4) 

where  p(x)  is  a  basis  (usually  polynomial)  and  a(x)  are  unknown  coefficients.  Some 
complete  polynomial  bases  are 


ID: 

•a 

>1 

II 

'  h-1' 

,v 

linear 

(5a) 

pT(x)  =  [l,x,*2] 

quadratic 

(5b) 

2D: 

II 

'h- 

linear 

(5c) 

pT(x)  =  [l,x,y,x2,xy,y2 

quadratic 

(5d) 

In  Section  5  it  is  shown  that  other  functions  can  be  added  to  the  basis  in  situations 
where  i);  is  desirable  to  enrich  the  solution. 

To  And  the  approximation  of  the  field  variable  by  Eq.  (4),  it  is  necessary  to  deter¬ 
mine  the  coefficients  a(x).  The  moving  least  squares  (MLS)  methodology  (Lancaster 
and  Salkauskas,  1981)  is  an  effective  technique  for  approximating  a  function  using  a 
set  of  scattered  data.  Given  a  set  of  nodes  with  coordinates  X/  at  which  the  field 
variable  u;  is  known,  a  weighted  £2  norm  can  be  written 

n 

J  =  tt»(x  -  X/)  [pT(x/)a(x)  -  U/]2  .  (6) 

/= 1 

where  w(x  —  X/)  is  the  weight  function  of  node  I  at  point  x,  and  n  is  the  number  of 
neighbors  of  point  x,  i.e.,  nodes  with  iv(x  —  X/)  >  0. 

Finding  the  minimum  of  J  with  respect  to  a(x)  leads  to  a  set  of  linear  equations 

A(x)a(x)  =  C(x)u  (7) 
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where 


A(x)  =  ^  w(*  ~  x/) p(x/)pT(x/)  (8a) 

/=1 

C(x)  =  [u)(x-x1)p(xi),w(x-x2)p(x2),...  , w(x  —  x„)p(x„)]  (8b) 

u  =  [ui,u2,...  ,un] .  (8c) 

The  matrix  A(x)  is  often  called  the  moment  matrix  and  C(x)  is  a  matrix  of  column 
vectors.  Eq.  (7)  can  be  solved  for  a(x)  to  yield 

a(x)  =  A-1(x)C(x)u 

which  can  be  substituted  into  Eq.  (4)  to  yield  an  approximation  in  terms  of  the  nodal 
coefficients 


uh(x)  =  Y2  PT(x)A  1(x)C/(x)u;,  (9) 

7=1 

where  C/(x)  is  the  Ith.  column  of  C(x)  and  uj  is  the  nodal  coefficient  for  the  Ith 
neighbor  of  x.  Defining  the  shape  function  </>/(x)  as 

<MX)  =  pT(x)  A-1  (x)C/(x)  (10) 

allows  the  approximation  to  be  written  as 

/ 

7  n 

uh(x)  =  Y2^^ui'  (n) 

i- 1 

which  is  a  form  familiar  to  those  with  a  finite  element  background. 

The  spatial  derivatives  of  the  shape  functions,  computed  by  the  chain  rule,  are 

<£/,.-(x)  =  p?(x)A_1(x)C/(x)  +  pr(x)  [A“1(x)C/(x)  +  A_1(x)C/ij(x)]  (12) 

where  A"1  =  —  A-1AjjA-1.  Note  that  the  second  term  in  Eq.  (12)  is  expensive  to 
compute  because  of  the  term  A”1 .  Nayroles,  Touzot,  and  Villon  (1992)  computed  only 
the  first  term  of  the  derivatives  which  results  in  the  inability  of  their  approximation 
to  satisfy  the  patch  test.  Krongauz  and  Belytschko  (1997a)  have  shown  that  the 
method  can  be  rendered  convergent  by  a  Petrov-Galerkin  formulation. 
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2.1  Fast  shape  function  and  derivative  computation 

The  number  of  operations  required  to  form  shape  functions  and  derivatives  can  be 
decreased  by  the  procedure  in  Belytscbko,  Krongauz,  Fleming,  Organ,  and  Liu  (1996) 
and  Fleming,  Chu,  Moran,  and  Belytscbko  (1997).  The  shape  function  in  Eq.  (10) 
can  be  written  as 

<f>i(x)  =  pT(x)A~1(x)Cj(x)  =  7T(x)C/(x)  (13) 

with  corresponding  derivatives 

<MX)  =  7j(x)c/(x)  +  7T(x)C/,,(x)  .  (14) 

Comparing  the  underlined  terms  in  Eq.  (13)  leads  to  the  relationship 

A(x)7(x)  =  p(x).  (15) 

The  coefficients  *y(x)  can  be  obtained  by  an  LU  decomposition  of  A(x)  and  backsub- 
stitution.  This  requires  fewer  computations  than  a  full  inversion  of  A(x),  which  is 
required  to  form  the  shape  functions  with  Eq.  (10). 

The  derivative  of  'y(x)  is  required  to  compute  the  shape  function  derivatives  in 
Eq.  (14).  Taking  the  derivative  of  Eq.  (15) 

-Mxb(x)  +  A(x)7it(x)  =  p,,-(x)  (16) 

and  rearranging  known  terms  which  to  the  right-hand  side  leads  to 

A(x)7,i(x)  =  P,.'(x)  ~  a,.(x)7(x)  •  (17) 

Using  the  LU  decomposition  of  A(x),  which  is  known  from  solving  Eq.  (15),  7,»(x) 
can  be  computed  with  only  a  backsubstitution.  Higher  order  derivatives  can  be  easily 
obtained  by  repeating  the  procedure  used  for  computing  the  first  derivative. 

While  this  procedure  for  computing  shape  functions  and  their  derivatives  is  theo¬ 
retically  identical  to  directly  evaluating  Eqs.  (10)  and  (12),  the  number  of  computa¬ 
tions  is  reduced.  In  addition  to  increasing  speed,  this  modification  can  also  alleviate 
ill  conditioning  of  shape  functions  in  cases  where  the  moment  matrix,  A(x),  is  not 
well  conditioned. 

2.2  Modeling  cracks  with  EFG 

A  crack  is  modeled  in  EFG  by  defining  a  line  segment  internal  to  the  domain.  The 
domains  of  influence  for  nodes  near  the  crack  are  truncated  whenever  they  intersect 
the  crack  surface  so  that  a  node  on  one  side  of  the  crack  will  not  affect  points  on  the 
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Crack 


New  crack  tip  node  — y* 
< 

\ 


Split  node 


\ 


Old  crack  tip  node 


Crack 


(a)  Initial  setup  (b)  After  1  step 

Figure  2:  Crack  setup  for  modeling  crack  propagation  with  EFG. 


opposite  side  of  the  crack.  This  technique  was  first  presented  in  Belytschko,  Lu,  and 
Gu  (1994)  has  been  called  the  visibility  criterion  by  Krysl  and  Belytschko  (1996). 
The  nodal  domain  of  influence  can  be"  considered  as  the  fine  of  sight  and  a  crack 
can  be  considered  as  an  opaque  boundary.  Whenever  the  line  of  sight  meets  the 
opaque  boundary,  the  domain  of  influence  is  cut.  The  visibility  criterion  has  some 
limitations  near  crack  tips  which  are  discussed  in  Section  4.  Determining  whether  a 
line  drawn  between  a  node  and  a  sampling  point  intersects  a  crack  segment  requires  an 
efficient  and  robust  algorithm  capable  of  handling  complex  geometries.  The  Sedgewick 
algorithm  (Sedgewick,  1988)  is  found  to  work  well. 

One  of  the  biggest  benefits  of  meshless  methods  such  as  EFG  is  their  inherent 
ability  to  model  arbitrary  crack  propagation  due  to  the  absence  of  a  predefined  ele¬ 
ment  connectivity.  Growing  a  crack  consists  of  simply  adding  another  line  segment 
at  the  existing  crack  tip  as  shown  in  Fig.  2. 

As  the  crack  propagates,  it  will  pass  directly  through  the  current  crack  tip  node 
and  may  pass  through  other  nodes  as  well.  In  this  case,  the  node  is  replaced  by  two 
nodes,  one  above  the  crack  and  another  below  the  crack  (see  Fig.  2).  This  is  preferred 
to  deleting  the  node  because  it  increases  the  spatial  resolution  along  the  crack  surface. 

Using  a  linear  basis  for  crack  problems  necessitates  increased  nodal  refinement 
around  the  crack  tip  in  order  to  capture  the  crack  tip  stress  field  accurately.  A  star¬ 
shaped  array  of  nodes  consisting  of  several  rings  of  nodes  around  the  crack  tip  is 
found  to  work  well.  A  recommended  arrangement  of  nodes  is  given  in  (Terry,  1994). 
The  star-shaped  arrangement  is  found  to  be  more  effective  if  existing  nodes  which 
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fall  within  the  radius  of  the  stax  are  turned  off  (except  for  the  crack  tip  node).  For 
a  propagating  crack,  the  star-shaped  array  moves  with  the  crack  tip  and  turned  off 
nodes  are  turned  back  on  once  they  are  no  longer  in  the  radius  of  the  star. 

3  Elastostatics 

The  report  focuses  on  fracture  of  linear  elastic  media.  The  elastostatic  boundary  value 
problem  is  reviewed,  the  variational  form  is  given  and  numerical  approximations  by 
mesh-free  methods  are  shown;  enforcement  of  essential  boundary  conditions  is  also 
discussed. 

Consider  a  two-dimensional  domain  fl  bounded  by  T.  The  equation  of  equilibrium 
is 

V  •  +  b  =  0  in  n  (18) 

where  &  is  the  stress  tensor  for  a  displacement  field  u,  and  b  is  the  body  force.  The 
boundary  conditions  are 

u  =  u  on  T„ 

<t  •  n  =  t  on  Tt 

where  the  superposed  bar  indicates  prescribed  values  and  n  is  the  unit  normal  vector 
to  T. 

The  variational  (or  weak)  form  for  Eq.  (18)  can  be  written 

5W±  f  Vs5u  :  <rdn  -  [  5n  ■  b dO  -  f  8u-tdT  -  <Wu(u)  =  0,  W5u  e  Ux  (19) 

Jn  Jn  Jrt 

where  Vs  is  the  symmetric  gradient  operator  and  5u  is  a  test  function  from  the  same 
space  of  functions  as  u,  i.e.,  this  is  a  Bubnov-Galerkin  method  (see  Hughes  (1987)). 
The  term  <5kFu(u)  is  required  for  enforcing  the  essential  boundary  conditions  in  a 
meshless  method  and  will  be  discussed  in  the  subsequent  section. 

For  lineax  elasticity, 

e  =  Vsu  <x  =  D  :  £  (20) 

which  can  be  used  to  write  the  weak  form  from  Eq.  (19)  in  terms  of  the  displacement, 

u. 

The  discrete  form  of  Eq.  (19)  for  a  meshless  method  can  be  obtained  using  the 
approximation  from  Eq.  (11)  as  an  approximation  for  u  and  <£u  to  get  the  system  of 
equations 


Ku  =  fext . 


(21) 
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The  stiffness  matrix  K  £  5ftne<?Xne<?  (neq  is  the  number  of  equation)  and  the  external 
force  vector  fext  £  3Rne<?  are  composed  of  the  submatrices 

K IJ=  [  B/DB jdSl  (22a) 

J  o 

ffxt  =  [  faidT  +  [  4>ibdfl  (22b) 

JTt  JO. 

where  K u  £  SRnsdXns<i,  fjxt  £  3?nsd  (nsd  is  number  of  spatial  dimensions),  D  is  the 
constitutive  matrix 

plane  strain  (23) 


plane  stress  (24) 


(25) 

3.1  ;Enforcement  of  essential  boundary  conditions 

f 

One  drawback  of  MLS  approximations  is  that  essential  boundary  conditions  cannot 
be  satisfied  directly  because  ^  Sjj,  and  consequently,  shape  functions  from 

nodes  on  the  interior  of  the  domain  are  nonzero  on  the  boundary.  The  term  51Tu(u) 
in  Eq.  (19)  is  one  way  of  enforcing  essential  boundary  conditions.  Some  forms  which 
have  been  suggested  are  1)  Lagrange  multipliers  (Belytschko,  Lu,  and  Gu,  1994) 

8WU{ u)  =  [  SX  ■  (u  -  u)dT  +  f  Su-XdT,  (26) 

J r„  J r„ 

where  A  is  a  Lagrange  multiplier,  2)  a  modified  variational  principle  in  which  the 
Lagrange  multipliers  are  replaced  by  their  physical  meaning,  the  traction  (Lu,  Be¬ 
lytschko,  and  Gu,  1994) 

JW„(u)  =  /  5t  •  (u  —  u)dT  +  f  <5u  •  tdT , 

J  r  ti  Jru 


D  = 


D  = 


E 


(1  +  I/)(l  —  2v) 
E 


l  —  i/  v  0 

v  l  —  i/  0 
0  0 


1  —  I/2 


1  1/  0 

1/  1  0 

0  0  ^ 


and  B/  is  a  matrix  of  shape  function  derivatives 


B  ,= 


4>i,x  o 
0  4>I,y 


(27) 
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where  t  =  <x  *  n,  and  3)  a  penalty  method  (Belytschko,  Gu,  and  Lu,  1994) 

2dT  (28) 

where  j3  is  a  penalty  parameter. 

Another  method  of  enforcing  essential  boundary  conditions  in  meshless  methods  is 
by  coupling  with  finite  elements  (Krongauz  and  Belytschko,  1996;  Belytschko,.  Organ, 
and  Krongauz,  1995).  In  this  method,  a  row  of  finite  elements  is  placed  along  the 
essential  boundaries  allowing  the  boundary  conditions  to  be  enforced  by  prescribing 
the  values  at  the  nodes. 

3.2  Integration  issues 

Computing  the  stiffness  matrix  and  force  vector  from  Eqs.  (22)  requires  an  integration 
over  the  domain  H,  which  in  two  dimensions  corresponds  to  an  area  integration.  Inte¬ 
grating  the  stiffness  matrix  and  force  vector  requires  a  numerical  integration  scheme 
such  as  Gauss  quadrature,  which  in  turn,  requires  a  subdivision  of  the  domain.  Since 
meshless  methods  have  no  inherent  subdivision  of  the  domain  like  finite  elements, 
it  is  necessary  to  introduce  a  subdivision  of  the  domain.  One  type  is  a  background 
cell  structure.  Two  subdivisions  are  shown  in  Figs.  3.  The  first  one  shown,  which 
is  the  most  common,  uses  a  finite  element  mesh  generator  to  create  a  cell  structure 
which  matches  the  domain;  this  technique  is  often  called  an  element  quadrature.  The 
vertices  of  this  background  mesh  are  often  used  as  the  initial  array  of  nodes  for  the 
EFG  model;  however,  additional  nodes  may  be  added  where  desired  such  as  the  nodes 
at  the  crack  tip  in  the  model  shown. 

The  second  integration  technique,  which  is  often  called  cell  quadrature ,  uses  a 
background  grid  of  cells  which  is  independent  of  the  domain.  At  each  integration 
point  it  is  necessary  to  determine  if  it  lies  inside  the  domain  before  it  is  used  for 
integrating  Eqs.  (22).  This  technique  is  not  widely  used  because  it  does  not  yield 
good  accuracy  along  curved  and  angled  boundaries. 

A  nodal  integration  technique  was  proposed  by  Beissel  and  Belytschko  (1996)  in 
an  effort  to  make  EFG  a  completely  meshless  method.  However,  additional  stability 
terms  must  be  added  to  make  the  method  stable  and  accuracy  is  not  as  good  as 
cell-based  integration  schemes. 

Hegen  (1997)  proposed  subdividing  the  cells  through  which  a  crack  passes  by  a 
triangulation  technique  to  avoid  integration  errors.  This  subdivision  is  not  found 
to  significantly  increase  the  accuracy  of  EFG  computations  and  imposes  undesirable 
“remeshing”  conditions  on  a  meshless  method.  The  lack  of  a  noticeable  effect  when 
aligning  the  cell  structure  with  a  crack  can  be  shown  by  using  the  model  shown  in 
Fig.  4.  A  crack  is  placed  between  rows  of  nodes  and  the  near  crack  tip  displacement 
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(b)  Cell  quadrature 

Figure  3:  Two  integration  methods  for  integrating  the  weak  form  with  a 
meshless  method. 
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Figure  4:  Discrete  model  for  near-tip  crack  problem. 


field  is  applied  to  the  boundary;  for  this  problem,  the  asymptotic  near-tip  field  is  the 
exact  solution.  Gauss  quadrature  is  performed  in  a  cell  structure  with  vertices  at 
the  nodes.  The  problem  is  also  solved  by  dividing  the  cells  through  which  the  crack 
passes  so  that  the  cell  boundaries  coincide  with  the  crack.  No  appreciable  change  in 
the  error  in  strain  energy  (1%  change)  or  stress  intensity  factor  (0.1%  change)  was 
noticed' when  the  crack  was  allowed  to  pass  through  the  middle  of  the  cell. 

In  this  report,  element  quadrature  is  used  with  4x4  Gauss  quadrature  points  in 
each  cell.  For  cells  near  a  crack  tip,  the  quadrature  order  is  increased  to  9  x  9.  Cells 
surrounding  a  crack  are  not  subdivided  to  align  cell  boundaries  with  the  crack. 

3.3  Domains  of  influence 

Properly  choosing  the  domain  of  influence  or  nodal  support  is  an  important  aspect 
of  meshless  methods.  The  size  of  the  support  should  be  sufficiently  large  so  that  the 
moment  matrix  is  regular  and  well  conditioned  and  that  the  spatial  distribution  of 
neighbors  is  fairly  even.  On  the  other  hand,  choosing  domains  of  influence  which  are 
too  large  leads  to  a  great  deal  of  computational  expense  in  forming  the  approximations 
as  well  as  assembling  the  stiffness  matrix.  Support  sizes  which  are  too  large  also 
detract  from  the  local  character  of  the  approximation;  for  problems  involving  sharp 
gradients,  some  loss  of  accuracy  is  typically  noted  as  the  effect  of  the  gradient  is 
smeared. 


4  SMOOTHING  OF  EFG  APPROXIMATIONS 


15 


These  aspects  of  meshless  methods  are  easily  shown  by  considering  a  one  dimen¬ 
sional  bar  on  the  domain  0  <  x  <  n  with  Young’s  modulus  of  E  =  1.  The  bar  is 
loaded  with  a  body  force  b(x)  =  sin(x)  and  boundary  conditions  u(0)  =  0,  =  -1. 

Evenly  spaced  nodal  arrangements  ranging  from  21  to  201  nodes  are  used  and  the 
characteristic  nodal  spacing,  c/,  is  set  equal  to  the  nodal  spacing;  the  nodal  support 
size  is  dmi  =  dmax  cj.  Fig.  5  presents  the  stress  profile  along  the  length  of  the  bar 
computed  using  dmax  =  2.0  and  dmax  =  1.01  and  Fig.  6  shows  results  for  the  er¬ 
ror  in  energy  computed  with  varying  support  sizes.  When  the  support  size  is  only 
slightly  larger  than  the  nodal  spacing,  the  error  is  high  and  the  stress  field  mimics 
a  linear  finite  element  solution,  i.e.,  the  stress  is  nearly  discontinuous  at  the  nodes. 
Increasing  the  support  size  to  dmax  =  2.0  or  2.5  leads  to  a  sharp  improvement  in 
accuracy,  although  the  rate  of  convergence  remains  the  same.  Further  increasing  to 
dmax  =  3.0  and  3.5  actually  leads  to  an  increase  in  error,  most  likely  because  the  gra¬ 
dients  in  the  solution  become  smeared  for  such  large  nodal  supports.  Note  that  the 
rate  of  convergence  in  each  case  is  linear  and  only  the  absolute  accuracy  increases  with 
increasing  dmax.  These  results  disagree  with  the  results  of  Liu,  Li,  and  Belytschko 
(1996)  and  Duarte  and  Oden  (1996a),  who  showed  increasing  rates  of  convergence. 
This  increasing  rate  of  convergence  is  perhaps  attributable  to  keeping  the  domain 
of  influence  constant  in  size  as  the  mesh  is  refined.  Note  that  the  RKPM  in  Liu, 
Li,  and  Belytschko  (1996)  and  the  hp~ clouds  both  used  the  moving  least  squares 
approximations  (see  Belytschko,  Krongauz,  Organ,  Fleming,  and  Krysl  (1996)). 

Determining  proper  domains  of  influence  in  two  and  three  dimensions  is  more 
difficult  than  in  one  dimension.  The  stress  fields  are  usually  much  more  complicated 
and  evenly  spaced  nodal  arrangements  are  not  practical. 

f 

4  Smoothing  of  EFG  approximations  near  noncon- 
vex  boundaries 

The  smoothness  which  is  inherent  in  meshless  methods  is  a  two-edged  sword.  On 
one  hand,  it  provides  approximations  of  functions  and  their  derivatives  which  are 
smooth  and  have  the  same  continuity  as  the  weight  function.  However,  in  cases 
where  a  discontinuity  is  present  in  either  the  geometry  or  the  material,  this  higher 
order  smoothness  leads  to  difficulties  which  must  be  addressed  in  order  to  obtain 
good  accuracy.  Discontinuities  can  arise  from  material  and  geometric  sources.  A 
material  discontinuity  occurs  along  an  interface  between  two  materials  with  differ¬ 
ent  properties,  leading  to  discontinuous  strains  across  the  interface.  This  situation 
is  easily  modeled  using  finite  elements  where  the  interpolants  are  C°.  Cordes  and 
Moran  (1996)  have  solved  problems  with  multiple  materials  using  EFG  by  treating 
the  individual  materials  as  separate  bodies  and  joining  them  together  with  Lagrange 
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x  postition 

Figure  5:  Stress  distribution  along  one  dimensional  bar  using  20  nodes. 


Figure  6:  Error  in  energy  versus  nodal  spacing  with  various  support  sizes  for 
a  one-dimensional  bar  with  b(x)  =  sin(x). 
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multipliers.  Krongauz  and  Belytschko  (1997b)  have  developed  techniques  in  which  a 
so-called  jump  term  is  included  in  the  trial  function  which  is  capable  of  representing 
the  discontinuity.  The  magnitude  of  the  jump  becomes  an  unknown  in  the  Galerkin 
discretization. 

The  second  class  of  problems  which  contain  discontinuities  are  due  to  geometrical 
effects.  This  includes  cases  where  a  boundary  of  the  geometry  is  nonconvex,  such 
as  a  body  with  a  hole  or  a  crack.  Because  of  the  aforementioned  smoothness  of 
meshless  methods,  special  procedures  are  required  near  nonconvex  boundaries.  In 
this  section,  techniques  for  modeling  nonconvex  boundaries  will  be  discussed.  The 
visibility  criterion  will  be  presented  first  and  some  of  its  limitations  will  be  pointed 
out.  To  overcome  some  of  the  limitations  of  the  visibility  criterion,  the  diffraction , 
transparency  and  “see-through”  methods  are  introduced.  Numerical  examples  will  be 
given  to  illustrate  the  behavior  of  these  alternative  techniques  and  show  when  they 
are  necessary. 


4.1  Visibility  Criterion 

The  first  technique  for  dealing  with  nonconvex  boundaries  is  the  visibility  criterion 
(Belytschko,  Lu,  and  Gu  (1994),  the  name  was  coined  in  Krysl  and  Belytschko  (1996)). 
In  this  straightforward  approach,  the  domain  of  influence  is  considered  as  the  field  of 
vision  at  a  node.  All  boundaries,  internal  and  external,  are  considered  to  be  opaque 
so  that  the  field  of  vision  is  interrupted  as  soon  as  such  a  boundary  is  encountered. 
Consider  node  J  in  Fig.  7,  where  the  surface  of  the  crack  is  within  its  domain  of 
influence  and  is  therefore  truncated.  This  truncation  will  create  a  discontinuity  in 
the  shape  function  for  node  J  which  will  lead  to  the  desired  discontinuity  in  the 
solution  across  the  crack  (see  Fig.  8a, c). 

A  difficulty  with  the  visibility  criterion  arises  for  nodes  near  the  end  of  a  discon¬ 
tinuity,  i.e.,  a  crack  tip  in  fracture  mechanics.  Consider  node  I  in  Fig.  7.  The  field  of 
vision  is  cut  by  the  crack,  leading  to  a  discontinuity  along  line  AC,  i.e.,  the  line  of  the 
crack.  However,  the  field  of  vision  is  also  truncated  along  line  AB,  which  extends  into 
the  domain.  This  leads  to  a  discontinuity  in  the  weight  function  as  well  as  the  shape 
function  along  this  line  as  shown  in  Fig.  8b, d.  Note  that  because  the  shape  functions 
are  created  from  the  interaction  of  weight  functions,  additional  discontinuities  arise 
in  the  shape  functions  from  other  nodes  which  have  discontinuous  weight  functions 
due  to  the  crack  tip. 

Discontinuities  in  the  shape  functions  can  also  arise  for  nodes  near  the  surface  of 
nonconvex  boundaries  such  as  a  hole.  As  an  example,  consider  the  domain  of  influence 
of  node  I  in  Fig.  9.  Whenever  a  ray  emanating  from  the  node  grazes  the  surface  of  the 
hole  such  as  points  A  or  C,  the  weight  function  will  have  a  discontinuity  along  lines 
AB  and  CD.  A  contour  plot  of  the  weight  function  for  node  I  is  shown  in  Fig.  10a 
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DOMAIN  OF  INFLUENCE  DOMAIN  OF  INFLUENCE 

FOR  NODE  J  FOR  NODE  I 


Figure  7:  Domain  of  influence  by  the  visibility  criterion  near  a  crack. 


and  the  corresponding  shape  function  is  shown  in  Fig.  10b.  Note  that  the  shape 
function  has  several  lines  of  discontinuity  because  it  is  formed  by  the  interaction  of 
several  weight  functions,  some  of  which  have  discontinuities  similar  to  those  shown 
in  Fig.  10a. 

The  presence  of  discontinuities  within  the  domain  is  undesirable  in  a  Galerkin 
method  and  must  be  handled  with  care.  The  length  and  size  of  the  discontinuities  is  a 
function  of  the  nodal  refinement  near  a  nonconvex  boundary,  i.e.,  as  the  nodal  spacing 
goes  to  zero,  the  discontinuities  go  to  zero.  Using  this  argument  and  others  similar  to 
the  arguments  used  in  proving  convergence  for  nonconforming  finite  elements,  Krysl 
and  Belytschko  (1996)  showed  that  the  discontinuous  approximations  generated  by 
the  visibility  criterion  can  lead  to  convergent  solutions. 


4.2  Diffraction  Method 

Continuous  and  smooth  approximations  near  nonconvex  boundaries  can  be  con¬ 
structed  quite  easily  by  the  diffraction  method  (Organ,  Fleming,  and  Belytschko, 
1996;  Organ,  1996;  Belytschko,  Krongauz,  Fleming,  Organ,  and  Liu,  1996).  The 
nodal  support  is  wrapped  around  nonconvex  boundaries  similar  to  the  way  light 
diffracts  around  sharp  corners.  This  method,  which  has  also  been  called  the  wrap¬ 
around  method,  is  quite  general  and  can  be  used  for  cracks  or  smooth  boundaries  such 
as  interior  holes.  The  method  will  first  be  described  for  cracks  and  then  generalized 
for  other  nonconvex  boundaries. 
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(c)  Shape  function  for  node  J  (d)  Shape  function  for  node  I 

Figure  8:  Contour  plots  of  the  weight  function  tn(x  —  X/)  and  shape  function 
as  determined  by  the  visibility  criterion  for  nodes  adjacent 
to  a  line  of  discontinuity  due  to  a  crack. 
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Figure  9:  Domain  of  influence  near  an  interior  hole  by  the  visibility  criterion. 


(a)  Weight  function  for  node  A  (b)  Shape  function  for  node  A 

Figure  10:  Contours  for  weight  and  shape  function  near  an  interior  hole  by 
the  visibility  criterion. 
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(a)  Support  around  a  crack  tip. 


Figure  11:  The  diffraction  (wrap-around)  method  for  constructing  smooth 
weight  functions  around  nonconvex  boundaries. 


Consider  Fig.  11a,  where  a  line  between  the  node,  X/,  and  a  sampling  point,  x, 
intersects  a  crack  and  the  tip  is  within  the  domain  of  influence  of  the  node.  The 
weight  function  distance,  dj ,  is  modified  (lengthened)  by 

(29) 

where  sx  =  ||xj  -  xc||,  s2(x)  =  ||x  -  xc||,  s0(x)  =  ||x  -  x/||,  and  xj  is  the  node,  x  is 
the  sampling  point,  and  xc  is  the  crack  tip.  The  parameter  A  is  used  to  adjust  the 
distance  of  the  support  on  the  opposite  side  of  the  crack.  It  was  found  that  A  =  1, 2 
performs  well.  Contour  plots  of  the  weight  and  shape  functions  by  the  diffraction 
method  are  shown  in  Fig.  12  and  surface  plots  are  shown  in  Fig.  13. 

The  spatial  derivatives  of  the  weight  function  are  computed  using  the  chain  rule: 

dw  _  dw  ddi 
dxi  ddj  dxi 

Since  dw/ddj  is  unchanged,  all  that  is  necessary  are  expressions  for  ddj/dxi : 

p  =  \  (S-L±n)x"  pL  +  (1  _  A)  (iL+iiV  pt, 
dxi  \  s0  J  ox;  \  s0  J  dxi 


(31) 
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(c)  Weight  function  for  A  =  2.  (d)  Shape  function  for  A  =  2. 

Figure  12:  Contours  for  weight  and  shape  functions  associated  with  node  A 
near  a  crack  tip  constructed  using  the  diffraction  method.  The 
quartic  weight  function  in  (lb)  is  used  with  dmax  =  2.01. 


4.2  Diffraction  Method 


(a)  Weight  function  associated  with  node 
A  in  Fig.  12a. 


z 


(b)  Shape  function  associated  with  node 
A  in  Fig.  12b. 

Figure  13:  Surface  plots  for  weight  and  shape  functions  near  a  crack  tip  con¬ 
structed  using  the  diffraction  method  with  A  =  1.  The  quartic 
weight  function  in  (lb)  is  used  with  dmax  —  2.01. 
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where 


dso 

dxi  s0 


and 


9s2  '  Xi 
dxi  s2 


The  diffraction  method  works  well  for  general  nonconvex  boundaries  as  well.  Con¬ 
sider  the  case  shown  in  Fig.  lib  where  the  line  between  a  node  and  a  sampling  point 
intersects  the  boundary  of  a  hole.  The  tangent  point  between  the  node  and  the  non¬ 
convex  boundary  is  used  as  the  wrap-around  point,  xc,  and  Eq.  (29)  used  to  compute 
the  weight  function  distance,  d/  (see  Fig.  lib).  Note  that  for  this  implementation 
of  the  diffraction  method,  the  segment  d2(x)  intersects  the  boundary.  This  should 
not  be  a  cause  for  alarm  because  with  adequate  refinement  this  effect  is  quite  negligi¬ 
ble.  As  is  shown  in  Section  4.4,  acceptable  results  can  be  obtained  when  the  smooth 
boundary  is  considered  fully  transparent. 


4.3  Transparency  Method 

Another  technique  for  constructing  continuous  approximations  is  the  transparency 
method  (Organ,  Fleming,  and  Belytschko,  1996;  Belytschko,  Krongauz,  Fleming,  Or¬ 
gan,  and  Liu,  1996),  which  will  be  described  here  for  cracks.  The  underlying  concept 
of  this  method  is  to  endow  the  crack  tip  with  a  varying  measure  of  transparency  such 
that  it  is  completely  transparent  at  the  tip  and  becomes  completely  opaque  a  short 
distance  behind  the  tip.  In  this  way,  the  field  of  vision  for  a  node  near  the  crack  tip  is 
not  abruptly  truncated  when  it  reaches  the  crack  tip,  but  rather  diminishes  smoothly 
to  zero  a  short  distance  behind  the  tip  of  the  crack. 

When  a  ray  passes  between  a  node  X/  and  a  sampling  point  x,  and  crosses  the 
crack  as  shown  in  Fig.  14,  the  distance  parameter  dj  in  the  weight  function  is  modified 
(lengthened)  by  the  following: 

d/(x)  =  s0(x)  +  dmI  ^  >  2,  (32) 

where  so(x)  =  ||x  =  X/||,  dmj  is  the  radius  of  support  for  node  I,  and  sc(x)  is  the 
intersection  distance  behind  the  crack  tip.  The  parameter  sc  sets  the  distance  behind 
the  crack  tip  at  which  complete  opacity  occurs: 

sc  =  nh ,  (33) 

where  h  is  the  nodal  spacing  and  k  is  a  constant,  usually  0  <  k  <  1. 

The  spatial  derivatives  of  the  distance  parameter,  d; .  obtained  by  the  chain  rule, 
are 

ddi  _  ds0  5*-1 

dxi  dxi  mI  s *  dxi 


(34) 
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Figure  14:  Transparency  method  for  computing  smooth  weight  functions. 


where 


Osq  _  XI{ 

dxi  so 


dsc  Xb  xc 

— —  =  —  cos  v  = - 

OX  1  sc 


and 


dsg 

dx2 


—  sin# 


Vb  -  yc 

$c 


where  Q  is  the  angle  between  the  crack  and  the  x-axis  and  x*,  is  the  intersection  point 
behind  the  crack  tip. 

Contour  plots  of  the  weight  and  shape  functions  near  a  crack  tip  constructed  by 
the  transparency  method  are  shown  in  Fig.  15  and  surface  plots  are  shown  in  Fig.  16. 
Note  that  the  functions  axe  continuous  at  the  crack  tip. 

One  drawback  of  the  transparency  method  is  that  it  does  not  work  well  when 
nodes  are  placed  too  close  to  the  crack  surface.  Fig.  17  shows  a  surface  plot  of  a 
shape  function  constructed  by  the  transparency  method  when  nodes  are  placed  along 
the  crack  surface.  Note  the  trough  which  appears  in  the  shape  function  ahead  of  the 
crack.  This  trough  appears  because  although  the  crack  tip  is  transparent  for  this 
node,  the  change  in  the  degree  of  transparency  with  respect  to  the  change  in  angle  is 
very  shaxp,  i.e.,  sc(x)  in  Eq.  (32)  increases  rapidly.  There  is  no  discontinuity  in  the 
shape  function,  only  a  small  dip  in  the  shape  function.  To  circumvent  this  problem  in 
the  transparency  method,  a  restriction  has  been  placed  on  the  position  of  the  nodes. 
All  nodes  should  be  placed  so  that  the  normal  distance  from  the  node  to  the  crack 
surface  is  greater  than  roughly  1/4 h,  where  h  is  the  nodal  spacing. 
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(c)  Weight  function  for  k  =  0.5.  (d)  Shape  function  for  k  =  0.5. 

Figure  15:  Contours  for  weight  and  shape  functions  associated  with  node  A 
near  a  crack  tip  constructed  using  the  transparency  method.  The 
quartic  weight  function  in  (lb)  is  used  with  dmax  =  2.01. 
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(a)  Weight  function  associated  with  node 
A  in  Fig.  15a. 
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(b)  Shape  function  associated  with  node 
A  in  Fig.  15b. 


Figure  16:  Surface  plots  for  weight  and  shape  functions  near  a  crack  tip  us¬ 
ing  the  transparency  method  with  k  =  1.0.  The  quartic  weight 
function  in  (lb)  is  used  with  dmax  =  2.01. 
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Figure  17:  Surface  plot  for  shape  function  using  the  transparency  method 
(k  =  1)  when  nodes  are  placed  too  close  to  the  crack  surface. 


4.4  “See-through”  Method 

Terry  (1994)  constructed  a  smooth  approximation  near  nonconvex  boundaries  using  a 
11  see-through”  method  in  which  all  or  part  of  the  curved  boundary  is  made  completely 
transparent.  Terry  (1994)  found  that  better  accuracy  is  obtained  for  a  problem  with 
an  interior  hole  when  the  boundary  of  the  hole  was  not  strictly  enforced  by  the 
visibility  criterion. 

Duarte  and  Oden  (1996b)  and  (Krysl  and  Belytschko,  1996)  suggested  a  smoothing 
technique  in  which  the  crack  was  completely  transparent  if  the  crack  tip  is  within  the 
domain  of  influence  of  a  node.  This  is  also  called  the  continuous  line  criterion :  if 
a  continuous  line  connecting  the  node  to  a  point  lies  entirely  withing  the  domain  of 
influence  of  the  node,  the  point  is  visible.  While  this  technique  is  easy  to  implement 
and  provides  smooth  approximations,  it  effectively  shortens  the  crack  and  leads  to 
inaccurate  solutions  (see  Fig.  18).  One  situation  in  which  this  method  does  work 
for  cracks  is  when  the  enrichment  techniques  from  Chapter  5  are  used.  In  this  case, 
the  approximation  function  near  the  crack  tip  is  strong  enough  to  yield  correct  crack 
opening  displacements. 
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Figure  18:  Crack  opening  displacement  and  Mises  stress  contours  when  us¬ 
ing  the  “ see-through ”  method  ( dmax  =  2).  The  dots  are  the  nodal 
locations  and  the  contours  are  generated  using  values  at  the  inte¬ 
gration  points. 


4-5  Numerical  Examples 

Two  problems  are  presented  here  to  illustrate  the  performance  of  the  smoothing 
techniques  for  nonconvex  boundaries.  The  EFG  method  is  used  for  numerical  com¬ 
putations  and  a  background  element  mesh  is  used  for  integrating  the  weak  form. 

4.5.1  Infinite  plate  with  a  hole 

An  infinite  plate  with  a  hole  subjected  to  a  remote  unit  traction  in  the  x-direction  is 
solved  and  the  solutions  are  compared  for  the  visibility,  diffraction,  and  “see-through” 
methods.  The  solution  of  this  problem  is  given  in  Timoshenko  and  Goodier  (1970) 
as: 


crxx-l-  \  (|  cos(20)  +  cos  (40))  +  ^  cos(40), 
ayy  =  cos(20)  -  cos(40))  -  cos (40), 

n2  1  8o4 

o-xy  =  “(^  sin(20)  +  sin(40))  +  —  sin(40), 


where  a  is  the  radius  of  the  hole. 


(35a) 

(35b) 

(35c) 
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Figure. 19:  Typical  mesh  for  the  infinite  plate  with  a  hole  problem.  Due  to 
symmetry,  only  a  quarter  mesh  is  modeled. 


A  finite  model  of  the  problem  is  constructed  by  applying  the  exact  tractions 
corresponding  to  Eq.  (35)  on  the  boundaries;  due  to  symmetry,  only  a  quarter  of  the 
plate  is  modeled  (see  Fig.  19).  The  dimensions  used  are  a  =  1.0  in.,  /  =  5.0  in.;  plane 
strain/ linear  elastic  conditions  are  assumed  with  Young’s  modulus  and  Poisson’s 
ratio,  E  =  30  x  106  psi  and  i /  =  0.3,  respectively.  A  typical  nodal  arrangement  is 
shown  in  Fig.  19.  The  EFG  meshes  were  graded  with  additional  refinement  around 
the  hole  up  to  a  radius  of  2a;  a  typical  mesh  is  shown  in  Fig.  19.  The  integration 
cells  coincided  with  the  nodal  arrangement;  9x9  Gauss  quadrature  is  used  in  the 
graded  region,  with  5x5  quadrature  in  the  remaining  cells;  for  the  simulations  by 
the  diffraction  method,  A  =  2. 

According  to  the  analytic  solution  in  Eq.  (35a),  the  hole  acts  as  a  stress  raiser 
with  a  stress  concentration  factor  of  3  at  ( x,y )  =  (0,a).  The  numerical  solution 
captures  the  stress  concentration  well  for  smaller  support  sizes  (see  Fig.  20),  but  the 
accuracy  deteriorates  for  larger  support  sizes  when  shape  functions  are  discontinuous. 
The  errors  in  the  stress  concentration  using  the  visibility  criterion  at  large  dmax  can 
be  traced  to  oscillations  in  the  stress  distribution,  as  shown  along  the  line  x  =  0  in 
Fig.  21. 

The  error  in  energy  as  a  function  of  the  domain  of  influence  is  plotted  in  Fig.  22, 
where  the  exact  solution  is  computed  from  the  solution  given  in  Timoshenko  and 
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Figure  20:  Convergence  of  the  stress  concentration  at  the  hole  versus  the 
number  of  nodes. 


Goodier  (1970),  and  the  energy  norm  is  computed  by 

1/2 

(36) 

The  results  again  show  that  for  smaller  domains  of  influence,  there  are  no  significant 
differences  in  accuracy  between  the  discontinuous  and  the  smooth  approximation 
functions;  however,  as  the  domain  of  influence  is  increased,  the  error  increases  for  the 
discontinuous  shape  functions.  On  the  other  hand,  both  accuracy  and  convergence 
rate  increase  for  the  smooth  shape  functions  as  dmax  increases  (see  Fig.  23). 

4.5.2  Mode  1  static  fracture 

A  near-tip  crack  problem  is  studied  for  the  effects  of  smoothing  the  approximation 
for  crack  problems  and  shows  the  improvements  in  accuracy  when  using  the  EFG 
enrichment  techniques  from  Section  5. 

For  the  numerical  model,  the  asymptotic  mode  1  displacement  field  in  Eqs.  (60)  is 
applied  to  the  boundaries  of  an  EFG  mesh,  with  the  stress  intensity  factor  prescribed 
as  ki  =  1  psi-v/m.  The  mesh  dimensions  are  2a  x  2a,  with  a  crack  of  length  a  oriented 
such  that  the  crack  tip  was  located  at  the  center  as  in  Fig.  24.  Numerical  results  are 


energy  norm  =  —  f  (cr  —  crh)  :  (e  —  eh)dtt 
_2  Jn 
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Figure  22:  Error  in  energy  versus  the  support  size  dmax  for  the  infinite  plate 
with  a  hole  problem. 


Number  of  nodes 


Figure  23:  Convergence  of  the  error  in  energy  versus  the  number  of  nodes 
for  the  infinite  plate  with  a  hole  problem. 
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Figure  24:  Typical  mesh  for  the  mode  1  fracture  problem. 


presented  for  both  a  linear  basis  and  an  intrinsically  enriched  basis  from  Eq.  (50).  In 
addition,  the  effects  of  local  crack  tip  refinement  are  investigated. 

Fig.  25a  shows  the  circumferential  distribution  of  the  displacement;  note  the  dis¬ 
continuities  near  the  crack  tip  for  shape  functions  based  on  the  visibility  criterion. 
These  discontinuities  occur  at  fixed  angles  and  vanish  away  from  the  crack  tip.  In 
other  words,  the  discontinuities  occur  on  rays  from  the  nodes  to  the  crack  tip;  the 
weight  functions  associated  with  those  nodes  axe  discontinuous  along  these  rays.  Krysl 
and  Belytschko  (1996)  have  shown  that  the  solution  is  still  convergent  due  to  the  local 
nature  of  the  discontinuities. 

One  drawback  of  the  smooth  approximations  is  that  the  computed  crack  opening 
profile  is  not  parabolic  at  the  tip.  The  shape  functions  wrap  around  the  crack  tip 
and  the  crack  is  effectively  shortened  if  the  smoothing  effect  is  too  large  or  the  mesh 
is  too  coarse,  leading  to  the  crack  opening  displacement  (COD)  shown  in  Fig.  26a. 
The  maximum  stress  is  then  not  at  the  crack  tip,  but  is  shifted  a  small  distance  that 
depends  on  the  mesh  refinement.  This  effect  can  be  reduced  by  increasing  A  in  the 
diffraction  method,  or  decreasing  k  in  the  transparency  method.  When  the  mesh  is 
refined  locally  as  in  Fig.  26b  or  an  enriched  approximation  is  used  as  in  Fig.  26c,  the 
maximum  stress  shifts  to  the  crack  tip  and  the  COD  profile  becomes  parabolic. 

Smooth  approximations  provide  substantial  benefits  when  using  an  enriched  ap¬ 
proximation.  Fig.  27  shows  the  error  in  energy  in  (36)  plotted  as  a  function  of 
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(a)  Discontinuous,  visibility  criterion. 
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(b)  Continuous,  diffraction  method. 

Figure  25:  Angular  variation  of  displacement  around  the  crack  tip  at  various 
radial  distances.  The  jumps  at  45°,  90°,  and  135°  in  (a)  occur  due 
to  the  discontinuous  shape  functions  generated  by  the  visibility 
criterion. 


4.5  Numerical  Examples 


36 


(b)  Refinement  at  crack  tip  (linear  basis).  (c)  Uniform  mesh  (enriched  basis) 

Figure  26:  Crack  opening  displacement  and  Mises  stress  contours  for  a  mode 
1  fracture  problem  using  the  diffraction  method  (A  =  1).  The  dots 
are  nodal  locations  and  the  contours  are  generated  using  values 
at  the  integration  points. 
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Figure  27:  Convergence  of  the  error  in  energy  for  the  mode  1  crack  problem. 
For  the  diffraction  method,  A  =  2. 


Figure  28:  Convergence  of  the  mode  1  stress  intensity  factor  using  the  en¬ 
riched  basis.  For  the  diffraction  method,  A  =  2. 
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uniform  nodal  refinement  for  both  linear  and  enriched  bases  with  the  discontinuous 
and  smooth  approximations.  The  results  are  nearly  identical  for  continuous  and  dis¬ 
continuous  approximations  with  a  linear  basis  because  oscillations  are  the  dominant 
source  of  error.  The  enriched  basis  reduces  oscillations,  so  discontinuities  in  the  shape 
functions  become  a  noticeable  source  of  error.  Fig.  28  shows  the  convergence  of  the 
mode  1  stress  intensity  factor  using  the  enriched  basis.  The  stress  intensity  factor 
was  computed  using  the  J  integral  in  domain  form  (Moran  and  Shih,  1987).  For  the 
discontinuous  shape  functions,  the  stress  intensity  factor  is  low  by  as  much  as  6%  for 
a  coarse  mesh,  but  improves  as  the  mesh  is  refined.  For  the  smooth  approximations, 
the  stress  intensity  factor  is  within  0.1%  for  all  meshes. 

As  stated  previously,  in  spite  of  the  discontinuities  due  to  the  visibility  criterion, 
the  approximation  is  still  convergent.  Figs.  29  show  four  levels  of  increasing  nodal 
refinement  surrounding  the  crack  tip.  Fig.  30a  shows  the  error  in  energy  using  a  linear 
basis  plotted  as  a  function  of  the  number  of  nodes  added  to  the  crack  tip  region  and 
Fig.  30b  shows  the  corresponding  stress  intensity  factors  for  those  levels  of  refinement. 
It  is  readily  seen  that  as  under  these  circumstances  the  discontinuous  approximations 
do  not  degrade  the  accuracy  compared  to  the  smooth  approximations  for  either  total 
error  or  stress  intensity  factor.  In  addition,  the  smooth  approximations  generated  by 
the  “see-through”  method  are  greatly  in  error,  leading  to  the  conclusion  that  these 
approximations  should  not  be  used  in  conjunction  with  shaxp  nonconvex  boundaries 
such  as  cracks.  It  should  also  be  noted  that  over-refinement  of  the  crack  tip  region 
can  actually  lead  to  an  increase  in  error  if  the  mesh  away  from  the  crack  tip  is  not 
refined.  Fig.  30a  shows  that  for  the  highest  level  of  refinement  in  Fig.  29d,  the  error 
actually  increases  slightly  over  the  previous  refinement.  As  a  general  rule  for  meshless 
methods,  a  sharp  gradient  in  nodal  spacing  leads  to  error.  A  likely  source  of  this  error 
is  that  nodes  in  the  coarse  region  will  have  domains  of  influence  which  extend  into 
the  refined  region  while  the  converse  is  not  true. 

5  Enrichment  of  EFG  for  Crack  Tip  Fields 

Finite  element  methods  have  reached  a  high  degree  of  effectiveness  in  analysis  of  sta¬ 
tionary  cracks.  A  large  variety  of  approaches' have  evolved,  such  as  direct  application 
of  standard  elements,  singular  crack  tip  elements,  and  enriched  elements.  While  di¬ 
rect  application  of  standard  elements  is  the  most  straightforward,  a  high  degree  of 
mesh  refinement  is  required  near  the  crack  tip  to  capture  the  singular  stress  fields. 
Singular  crack  tip  elements  can  achieve  high  accuracy  in  linear  elastic  fracture  me¬ 
chanics  without  extremely  refined  meshes.  The  most  popular  of  these  elements  is 
the  quarter-point  element  (Henshell  and  Shaw,  1975;  Barsoum,  1977;  Banks-Sills  and 
Bortman,  1984)  which  introduced  a  singularity  in  the  Jacobian  of  an  8-node  serendip- 


5  ENRICHMENT  OF  EFG  FOR  CRACK  TIP  FIELDS 


39 


(a)  2  rings,  5  nodes/ring  (b)  4  rings,  9  nodes/ring 


(c)  6  rings,  13  nodes/ring  (d)  8  rings,  17  nodes/ring 

Figure  29:  Four  levels  of  nodal  refinement  using  a  star-shaped  nodal  array 
at  the  crack  tip. 
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(a)  Error  in  energy 


(b)  Stress  intensity  factors 

Figure  30:  Energy  error  and  stress  intensity  factors  with  four  crack  tip  nodal 
refinements.  All  calculations  were  made  with  a  linear  basis  and 
Gaussian  weight  function  ( dmax  —  2.5,  a  =  0.625). 
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ity  element  by  moving  tbe  mid-side  node  on  the  sides  connected  to  the  crack  tip  to 
the  quarter-point.  Akin  (1976)  presented  an  element  in  which  the  shape  function  at 
the  crack  tip  node  is  modified  so  it  possesses  singular  derivatives  of  the  correct  order. 

Another  method  of  capturing  the  singular  stress  fields  in  finite  elements  is  to  enrich 
isoparametric  finite  elements  by  including  the  near  tip  fields  in  the  trial  functions 
(Benzley,  1974;  Gifford  and  Hilton,  1978).  The  advantage  of  these  elements  is  that 
the  stress  intensity  factors  can  be  computed  directly  as  part  of  the  solution.  However, 
the  results  of  elements  enriched  by  singular  fields  axe  quite  sensitive  to  their  size,  and 
do  not  exhibit  uniform  convergence.  In  addition,  these  elements  axe  cumbersome 
to  implement  because  the  stiffness  matrix  and  force  vector  have  to  be  expanded  to 
account  for  the  extra  unknowns,  the  stress  intensity  factors.  Transition  elements  are 
needed  because  these  enriched  elements  are  generally  not  compatible  with  standard 
elements. 

In  applications  of  EFG  or  meshless  methods,  special  techniques  have  been  devel¬ 
oped  for  incorporating  the  singular  functions  associated  with  elastostatic  fracture  as 
an  alternative  to  large  arrays  of  nodes  at  the  crack  tip;  . the  latter  can  be  expensive 
and  awkward  for  problems  with  complex  geometry.  It  was  found  that  the  incorpo¬ 
ration  of  the  singular  fields  in  a  meshless  method  is  substantially  simpler  and  more 
trouble-free  than  in  finite  element  methods. 

Enrichment  of  a  meshless  method  may  be  carried  out  extrinsically  or  intrinsically. 
Extrinsic  enrichment  consists  of  adding  an  enrichment  function  to  the  trial  function. 
For  intrinsic  enrichment,  the  enrichment  functions  are  included  in  the  EFG  basis,  so 
that  no  additional  unknowns  are  needed  in  the  final  system  of  equations. 

In  this  section,  extrinsic  and  intrinsic  enrichment  techniques  are  presented  for 
meshldss  methods.  In  Sections  5.1  and  5.2,  two  methods  of  extrinsic  enrichment 
are  presented  and  discussed.  Section  5.3  presents  intrinsic  enrichment.  Techniques 
for  coupling  enriched  and  linear  approximations  are  discussed  in  Section  5.4  and  a 
mapping  for  kinked  cracks  is  given  in  Section  5.5.  Numerical  results  are  shown  in 
Section  5.6.  In  addition,  two  degrees  of  enrichment  are  discussed:  1)  full  enrichment 
by  all  linearly  independent  functions  of  the  asymptotic  near-tip  field,  and  2)  radial 
enrichment  consisting  of  only  the  y/r  field. 

5.1  Extrinsic  MLS  Enrichment 

In  extrinsic  enrichment  of  a  meshless  approximation,  a  function  closely  related  to 
the  solution  is  added  to  the  polynomial  expansion  of  MLS,  Eq.  (4).  For  example,  in 
linear  elastic  fracture  mechanics,  the  near  tip  asymptotic  field  or  its  constituents  can 
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be  added.  Tlie  approximation  takes  the  form 

uJ(x)  =  pT(x)att(x)  +  f;(J:;Oi<,(x)  +  ^(jLWl  «  =  1,2  (37) 

J=1 

where  it£(x)  denotes  the  approximation  for  ua(x),  p(x)  is  a  complete  polynomial 
basis  in  the  spatial  coordinates,  nc  is  the  number  of  cracks  in  the  model,  a<,(x)  are 
the  coefficients  of  the  polynomial  basis;  k{  and  k2  are  global  unknowns  associated 
with  crack  j.  Lower  case  Greek  subscripts  have  a  range  of  2  and  refer  to  Cartesian 
components.  Some  complete  bases  axe  given  in  Eqs.  (5). 

The  functions  Qlo(x)  and  Q2a(x),  which  describe  the  near-tip  displacement  field 
for  an  elastic  crack,  axe  (Williams,  1957;  Anderson,  1991) 

<3u<x) =^\f^cos  (I)  k  - 1 + 2sin2  G). 

0u(x) = G)  k + 1  _  w  G). 
«»w-5^-».G)  ;+i+2c°s2  G). 
q22(x) = -%^cos  G)  [k  ■ 1  ■ 2sin2  G). 

where  r  is  the  distance  from  the  crack  tip,  6  is  the  angle  from  the  tangent  to  the 
crack  path  at  the  crack  tip  (see  Fig.  31),  fJ.  is  the  shear  modulus  and  k  the  Kolosov 
constant  defined  as 

3  —  4i/  plane  strain 

(3  —  r/)/(l  +  v)  plane  stress. 

The  coefficients,  aa(x),  are  functions  of  the  spatial  coordinates  and  are  determined 
by  the  MLS  methodology  in  Section  2.  However,  additional  terms  arise  from  the 
inclusion  of  the  near-tip  field  and  so  the  MLS  formulation  will  be  rederived  here  in 
the  interest  of  completeness.  A  weighted,  discrete  norm  is  written 


where  n  is  the  number  of  points  in  the  neighborhood  of  x  for  which  the  weight 
function,  w(x  —  X/),  is  non-zero,  and  u;Q  is  the  a  component  of  the  nodal  value  at 


(38a) 

(38b) 

(38c) 

(38d) 
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Figure  31:  Local  coordinate  system  at  crack  tip. 


X/.  The  stationarity  of  J  with  respect  to  aa(x)  leads  to 

n  (  ne 

A(x)aa(x)  =  ^  ]  C/(x)  i  Uja  —  ^  ^2^2a(Xj)] 


(40) 


/=  1 


J  =  1 


where 

n 

A(x)  =  ^  io(x  -  x/)p(x/)pT(x/)  (41) 

i=i 

C/(x)  =  w(x  -  x/)p(x/) .  (42) 

It  should  be  noted  that  k{  and  k2  are  global  parameters  in  this  method  and  they  are 
considered  fixed  in  the  process  of  obtaining  the  parameters  aQ  for  the  local  fit. 
Solving  Eq.  (40)  for  a^x)  gives 

n  (  ne 

A_1(x)C/(x)  <  uIa  -  ^[k{Q 
i=i  l  j=i 

Expressing  Eq.  (37)  in  terms  of  the  nodal  parameter  uja  and  the  enriched  field  pa- 


a«(x)  = 


i«(x/)  +  klQLM] 


(43) 
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rameters  k[  and  kJ2  yields 

u«(x)  =  S<MX)  {  Ula  ~  £[^QL(x/)  +  HQLM] 


nc 


1=1 

nc 


i= 1 


(44a) 


+  5^<2i«(x) +  ^2«(x)l 


3=1 


1=1 

nc 


Ua(*)  =  ^<MX)“/« 

n 

Qi«(x)  “  5^/(x)0i«(x/) 

7=1 
n 

QL(x)  -  'E,-<h(*)QL(xi) 


+  E*i 
+  E*4 

i=i 


(44b) 


i=i 


where  the  shape  function,  <^>/(x),  is  defined  as 

<MX)  =  PT(x)A_1(x)C7(x).  (45) 

These  shape  functions  are  identical  to~  those  defined  in  Eq.  (10)  and  are  capable  of 
representing  the  smooth  part  of  the  solution.  Eq.  (44a)  will  be  written  as 


n  nc 

>£(*)  =  E  +  EWOLW  +  >4qL  (*)]>  (46) 

;  /=i  j=i 

where  the  modified  nodal  coefficients,  u/a,  are 

Ula  =  UJa  -  53[&iQia(Xj)  +  KQLMY  (47) 

3=1 

Note  that  the  enriched  finite  element  formulation  given  by  Benzley  (1974)  and 
Gifford  and  Hilton  (1978)  is  the  same  as  Eq.  (44b).  This  form  is  advantageous  in 
finite  elements,  but  not  necessary  in  a  meshless  method.  Using  the  approximation 
in  Eq.  (46)  is  simpler  to  implement  and  computationally  faster  because  the  terms 
S/=1  <A/(x)&iQia(x/)  and  ^"=1  4>i{x)k?2Q:2CS3ti )  are  not  subtracted  at  each  sampling 
point. 


5.2  Extrinsic  PU  enrichment 

Extrinsic  enrichment  of  meshless  methods  can  be  carried  out  using  partition  of  unity 
(PU)  methods  {Duarte  and  Oden,  1996b;  Melenk  and  Babuska,  1996;  Belytschko, 
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Krongauz,  Organ,  Fleming,  and  Krysl,  1996).  In  this  method,  the  approximation  is 
augmented  by  enrichment  functions  added  extrinsically  to  the  existing  EFG  approxi¬ 
mation  from  Eq.  (11).  This  basis  can  consist  of  higher  order  polynomials  or,  for  linear 
elastic  fracture  problems,  terms  from  the  asymptotic  near  tip  field  can  be  used.  The 
extrinsic  basis  is  smoothly  added  to  the  existing  approximation  by  multiplying  it  by 
a  partition  of  unity. 

The  essential  element  of  this  method  is  the  construction  of  a  partition  of  unity, 
which  can  be  obtained  by  the  MLS  methodology  (see  Section  2).  A  partition  of 
unity  <fi(x)h,  constructed  from  a  complete  polynomial  basis  of  order  k,  is  a  local 
approximation  for  which 

<£/(x)  =  1  •  (48) 

7=1 

It  can  easily  be  seen  that  MLS  approximations  are  partitions  of  unity  since  Eq.  (48) 
is  the  reproducing  condition  for  a  constant,  which  MLS  approximations  must  satisfy 
(see  Belytschko,  Krongauz,  Fleming,  Organ,  and  Liu  (1996)). 

Approximations  based  on  partitions  of  unity  take  the  form 

n  n  me 

uh(x)  =  ^/(x)u/ + Y2  52  (49) 

1=1  '  7=1  «=1 

where  u/  and  b!t  are  nodal  coefficients,  and  n  is  the  number  of  neighbors  of  point  x. 
The  vector  q(x)  is  called  the  extrinsic  basis  and  is  of  length  me.  For  enriching  linear 
elastic  fracture  problems,  this  basis  can  contain  the  y/r  radial  dependence  or  the  y/r 
radial  and  angular  6  dependence  which  spans  Eq.  (38).  A  superscript  k  is  added  to 
the  shape  functions  <pj  in  the  approximation  to  denote  the  polynomial  order  of  the 
basis  used  in  forming  the  partition  of  unity.  It  is  sometimes  convenient  to  construct 
the  partitions  of  unity  using  Shepard  functions  (i.e.,  k  =  1)  which  satisfy  constant 
consistency  and  add  the  terms  to  satisfy  linear  and  higher  order  consistency  to  the 
extrinsic  PU  basis. 

PU  methods  appear  to  provide  a  vehicle  for  local  enrichment.  The  partition 
of  unity,  <i>i(x)k,  can  be  formed  from  a  linear  basis  ( k  =  2),  which  yields  linear 
consistency.  Enrichment  of  the  approximation  may  be  carried  out  locally  by  adding 
the  known  form  of  the  solution  to  the  extrinsic  basis,  q(x),  at  nodes  in  the  region  in 
which  it  is  needed.  It  should  be  noted  that  the  enrichment  should  be  added  to  each 
node  whose  domain  of  influence  extends  into  the  region  to  be  enriched. 


5.3  Intrinsic  Basis  Enrichment 

Meshless  approximations  can  be  intrinsically  enriched  by  including  a  special  functions 
in  the  basis  (Flgming,  Chu,  Moran,  and  Belytschko,  1997).  For  example,  in  fracture 
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mechanics,  one  can  include  the  asymptotic  near-tip  displacement  field  in  Eqs.  (38),  or 
an  important  ingredient  such  as  y/r.  The  choice  of  functions  depends  on  the  coarse- 
mesh  accuracy  desired.  For  higher  accuracy,  include  the  full  asymptotic  field  from 
(38),  while  for  higher  speed  at  some  cost  of  accuracy,  only  the  y/r  function  can  be 
included  in  the  basis.  Both  methods  are  described  in  the  subsequent  sections. 

5.3.1  Full  enrichment 

In  full  intrinsic  enrichment  of  EFG  approximations  for  fracture  problems,  the  en¬ 
tire  neax-tip  asymptotic  displacement  field  is  included  in  the  basis.  Following  some 
trigonometric  manipulation,  it  can  be  shown  that  all  the  functions  in  Eqs.  (38)  are 
spanned  by  the  basis 

pT(x)  = 

(the  linear  terms  are  not  related  to  the  near-tip  fields  and  are  represented  through 
the  linear  completeness  of  the  EFG  approximant).  This  basis  can  be  used  in  Eq.  (4) 
and  leads  to  approximations  of  the  form 

uh(x)  =  pr(x)A~1(x)C/(x)  ui .  (51) 

/=1>  ' 

where  </>/(x)  is  the  enriched  EFG  shape  function. 

In  contrast  to  the  extrinsic  methods  presented  in  Sections  5.1  and  5.2,  this  method 
involves  no  additional  unknowns.  However,  because  of  the  increased  size  of  the  basis, 
additional  computational  effort  is  required  to  invert  the  moment  matrix,  A(x).  In 
addition,  the  domain  of  influence  must  be  enlarged  to  achieve  regularity  of  A(x).  For 
multiple  cracks,  four  additional  terms  would  have  to  be  added  to  the  basis  for  each 
crack;  a  method  for  coupling  an  enriched  basis  with  a  linear  basis  is  presented  in 
Section  5.4  which  avoids  this  difficulty. 

Using  an  enriched  basis  can  lead  to  a  moment  matrix  which  exhibits  some  ill- 
conditioning.  While  this  generally  does  not  affect  the  final  solution  when  the  enriched 
basis  is  used  globally,  it  does  have  a  deleterious  effect  when  trying  to  couple  shape 
functions  from  an  enriched  basis  to  shape  functions  from  another  basis  such  as  a 
linear  basis.  A  treatment  which  has  been  found  to  be  effective  for  dealing  with  ill- 
conditioning  is  to  reduce  the  number  of  computations  required  to  compute  the  shape 
functions  by  the  procedure  in  Section  2.1.  This  procedure  replaces  the  full  inversion  of 
the  moment  matrix  with  only  a  forward  reduction  and  back  substitution.  A  second 
method  found  .to  be  effective  for  dealing  with  ill-conditioning  is  to  diagonalize  the 


0  6  0  0 
1,  x,  y,  Vr  cos  \fr  sin  \fr  sin  -  sin  0,  Vr  cos  -  sin  8 


(50) 
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moment  matrix  by  Gram-Schmidt  orthogonalization  (Lu,  Belytsckko,  and  Gu,  1994), 
which  increases  the  linear  independence  of  the  moment,  equations. 

When  an  enriched  basis  is  used  at  any  node  of  a  mesh,  it  must  be  used  at  all  nodes 
or  a  special  technique  must  be  used  to  blend  it  to  nodes  with  a  different  basis.  Simply 
deleting  functions  from  the  basis  results  in  discontinuities  in  the  approximation. 

5.3.2  Radial  enrichment 

Meshless  approximations  can  be  partially  enriched  by  expanding  the  intrinsic  basis 
to  include  only  the  radial  variation  in  Eq.  (38).  This  leads  to  the  basis 

pr(x)  =  [i,x,y,  Vr]  (52) 

where  r  is  the  radial  distance  from  the  crack  tip.  This  partial  enrichment  is  feasible 
because  the  angular  variation  around  the  crack  tip  is  smooth,  but  the  radial  variation 
is  singular  in  the  stress. 

The  advantage  of  partial  enrichment  is  that  the  intrinsic  basis  is  only  expanded 
by  one  term  and  inverting  the  moment  matrix,  A(x),  to  form  the  shape  functions  is 
much  cheaper  than  for  full  enrichment.  In  addition,  it  does  not  seem  to  be  necessary 
to  use  the  smoothing  techniques  in  Chapter  4  because  there  are  no  discontinuities  in 
the  radial  direction;  the  discontinuities  in  the  angular  direction  lead  to  a  noticeable 
loss  of  accuracy  when  full  enrichment  is  used. 

5.4  Coupling  enriched  and  linear  approximations 

Enriching  the  approximation  for  the  entire  domain  of  a  problem  is  generally  unnec¬ 
essary  and  leads  to  unneeded  computational  expense.  For  example,  the  crack-tip 
singular  field  is  local  to  the  crack  tip  with  radius  ~  0.1a,  where  a  is  the  length  of  the 
crack.  Two  techniques  axe  presented  here  for  coupling  enriched  and  linear  approx¬ 
imations,  one  which  uses  a  consistent  coupling  to  maintain  C°  continuity  and  one 
which  does  not. 

The  first  technique  involves  coupling  the  approximation  over  a  transition  region 
as  a  linear  combination  of  the  enriched  linear  approximations  (this  is  similar  to  the 
way  Belytschko,  Organ,  and  Krongauz  (1995)  coupled  EFG  to  finite  elements).  The 
approximation  is  written 

ufc(x)  =  Ruenr(x)  +  (1  -  R)ulin(x).  (53) 

where  uenr(x)  is  the  enriched  approximation  and  ultn  is  the  linear  approximation;  R 
is  a  ramp  function  which  is  equal  to  unity  on  the  enriched  boundary  of  the  coupling 
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Linear 


Figure  32:  Schematic  for  the  coupling  of  enriched  and  linear  approximations. 


region  (r  =  ?'i )  and  equal  to  zero  on  the  linear  boundary  of  the  coupling  region 
(r  =  r2)  (see  Fig.  32).  Some  suggested  polynomial  ramp  functions  are 


{1  —  £  linear  ramp 

1  —  10£3  +  15£4  —  6£5  quintic  ramp 


(54) 


where  £  =  (r  —  ri)/(r2  —  rq),  r  is  the  radial  distance  from  the  crack  tip. 

The  coupled  approximation  for  the  intrinsically  enriched  basis  from  Section  5.3  is 
written 

uh(x.)  =  ^  <MX)U/  (55) 

1=1 


where 


Mx)  =  ofT(x)  +  (l  -  n)<f,‘r(x).  (56) 

and  <fijnr(x)  is  the  shape  function  formed  from  the  enriched  basis  of  Eq.  (50)  and 
(j>ljn{x)  is  the  shape  function  formed  from  a  linear  basis.  This  method  ensures  a  com¬ 
patible  displacement  field.  The  continuity  in  the  strain  field  depends  on  the  continuity 
of  the  ramp  function,  R  (i.e.,  the  linear  ramp  will  yield  continuous  displacements,  but 
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discontinuous  strains  at  r  =  ri,r2;  both  displacements  and  strains  will  be  continuous 
and  smooth  using  the  quintic  ramp). 

Through  numerical  experiments,  it  was  found  that  the  location  of  the  inner  radius 
of  the  ramp,  rq,  is  rather  arbitrary  and  can  actually  begin  at  the  crack  tip.  However, 
the  outer  radius,  r2,  needs  to  be  outside  the  singularity-dominated  zone  to  obtain 
good  accuracy. 

It  should  be  noted  that  it  is  possible  to  confine  the  enrichment  to  the  crack  tip 
region  by  simply  changing  from  an  enriched  to  a  linear  intrinsic  basis  away  from 
the  crack  tip.  The  distance  where  the  basis  is  changed  should  be  outside  the  stress 
intensity  factor  dominated  region,  with  an  outer  radius  similar  to  the  transition  region 
described  above.  One  drawback  of  this  method  is  that  the  displacement  field  at  the 
transition  point  will  be  C-1,  i.e.,  it  will  have  a  small  jump  in  displacement.  While  this 
is  theoretically  not  permissible  in  a  Galerkin  method,  it  is  acceptable  if  the  transition 
radius  is  sufficiently  far  from  the  crack  tip  so  that  jump  is  very  small.  In  this  case, 
the  small  jump  will  not  have  a  noticeable  effect  on  the  solution. 

These  two  methods  of  coupling  the  enriched  and  linear  approximations  were  found 
to  work  better  for  the  intrinsically  enriched  basis  in  Section  5.3  than  for  the  extrin¬ 
sic  MLS  enrichment  in  Section  5.1.  The  extrinsic  MLS  enrichment  is  sensitive  to 
discontinuities  and  some  loss  of  accuracy  was  noticed  for  this  coupling. 

Another  method  of  localizing  the  enrichment  to  the  crack  tip  region  is  to  use 
extrinsic  PUM  enrichment  from  Section  5.2.  The  extrinsic  basis  can  be  added  only 
to  nodes  near  the  crack  tip  region.  However,  it  is  necessary  to  enrich  all  nodes  which 
contribute  to  the  crack  tip  region  or  else  the  enrichment  is  incomplete  and  the  solution 
is  poor; 

5.5  Mapping  for  kinked  cracks 

The  enriched  field  of  Eq.  (38)  predicts  a  discontinuity  along  9  =  ±7r.  However,  for 
kinked  cracks,  this  discontinuity  falls  in  the  interior  of  the  domain  behind  the  crack 
tip.  Therefore,  a  mapping  is  required  to  align  this  discontinuity  in  enriched  field  with 
the  actual  crack  in  the  area  behind  the  leading  crack  segment  (i.e.,  the  area  to  the 
left  of  the  dotted  line  in  Fig.  33).  The  mapping  is  designed  so  that  a)  the  preceding 
segment  of  the  crack  is  rotated  to  be  aligned  with  the  leading  crack  segment,  b)  the 
field  along  the  line  perpendicular  to  the  leading  segment  (denoted  by  the  dotted  line) 
is  not  rotated,  and  c)  the  field  between  the  two  is  rotated  by  a  proportion  of  the  angle 
a  to  Or. 

We  define  an  angle  9r  as  the  angle  between  the  leading  crack  segment  and  the 
preceding  segment,  and  a  as  the  angle  between  the  leading  crack  segment  and  a 
line  from  the  kink  (or  rotation)  point  (xrot,  t/rot)  to  the  point  of  interest  (x,y).  Let 
9  =  a  —  9r  andx  =  y/(x  —  xrot)2  +  (y  —  yrot)2-  A  new  angle  9  is  defined 
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Figure  33:  Coordinate  mapping  used  for  kinked  cracks. 


where 


0  = 


0 

Ztt/2  —  Or 


0  7T 
0R  -  7t/2  2 


=  eca 


0  >  0 
0  <  0 


ca  = 


7r/2 

3tt/2  -  Tr 

7t/2 

Qr  —  n/2 


(57) 


which  i's  used  to  determine  the  coordinates  of  a  point  in  the  mapped  coordinate  system 


x*  =  xrot  —  r  cos  0  —  xtip  (58a) 

y*  =  2/rot  -  f  sin  0  -  ytip  ■  (58b) 

These  coordinates  are  used  to  compute  r*  and  0*  which  are  used  in  Eq.  (38)  or  Eq.  (50) 
to  compute  the  enriched  fields  at  (x,  y).  Note  that  the  computation  of  strain  involves 
derivatives  with  respect  to  (x*,y*)  and  it  is  necessary  to  obtain  the  derivatives  with 
respect  to  the  original  crack  tip  coordinate  system  (x,  y).  The  derivatives  are  obtained 
by  the  chain  rule  which  leads  to 


S1 

_ 1 

—  cos  0  cos  a  —  sin  0  sin  a  Ci  —  sin  0  cos  a  +  cos  0  sin  a  Ci 

d  I 

dx* 

1 _ 

—  cos  0  sin  a  +  sin  0  cos  aCi  —  sin  0  sin  a  —  cos  0  cos  a 

d 

.dy*. 

where  C/  =  Ca  x>r  Cb  depending  on  the  sign  of  0  (see  Eq.  (57)). 
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5.6  Numerical  Results  for  Enrichment 

Several  problems  are  solved  to  illustrate  the  effectiveness  of  enriching  the  EFG  formu¬ 
lation  for  fracture  problems.  Solutions  are  given  for  single  and  mixed  mode  problems. 

5.6.1  Near-tip  crack  problem 

A  closed  form  solution  for  a  crack  problem  can  be  constructed  by  using  the  well-known 
near-tip  field  in  a  domain  about  the  crack  tip  and  prescribing  the  displacements  along 
the  boundary  according  to  this  field.  This  is  called  a  patch  test  for  singular  fields. 
A  square  patch  with  sides  of  length  2 a  and  a  crack  of  length  a  is  used  (see  Fig.  34). 
This  problem  is  used  to  compare  the  performance  of  full  enrichment  (extrinsic  and 
intrinsic),  radial  enrichment,  and  a  linear  basis.  The  displacement  field  for  a  mode  1 
crack  is  (Williams,  1957) 

“*(x)  =  ^  Vs cos  (5)  [K ' 1 + 2  si“2  (5) .  (60a) 

“»<x)  =  sin  (5)  [x + 1 " 2cos2  (i)j  •  (60b) 

where  r  is  the  distance  from  the  crack  tip  and  6  is  the  angle  measured  from  the  line 
of  the  crack  (see  Fig.  31).  The  stress  intensity  factor  is  prescribed  as  k\  =  1  psi\/m. 
The  stresses  resulting  from  this  displacement  field  satisfy  equilibrium  so  the  solution 
is  exact  if  the  displacements  from  Eqs.  (60)  are  prescribed  on  the  outer  boundaries. 

Without  enrichment,  the  EFG  method  requires  considerable  nodal  refinement  near 
the  crack  tip  to  capture  the  singular  stress  field  with  sufficient  accuracy  to  provide 
a  useful  computation  of  the  stress  intensity  factor;  the  mesh  used  for  comparison  is 
shown  in  Fig.  34.  A  regular  grid  of  nodes  is  used  throughout  the  domain  with  a  radial 
array  of  nodes  around  the  crack  tip  for  enhanced  spatial  resolution.  The  background 
finite  element  mesh  used  for  quadrature  of  the  Galerkin  weak  form  is  shown  in  Fig.  34c. 
The  stress  field  ahead  of  the  crack  tip  computed  with  a  linear  basis  and  the  visibility 
criterion  is  shown  in  Fig.  35.  Without  the  radial  array  of  nodes  around  the  crack 
tip,  the  singular  stress  is  not  adequately  captured  by  regular  EFG  and  the  computed 
stress  intensity  factors  are  generally  5  —  15%  low  (see  Fig.  42).  The  radial  array  aids 
in  capturing  the  singular  stress  field;  however,  the  stresses  are  oscillatory  in  the  radial 
direction  (see  Fig.  35)  and  in  the  angular  direction  (see  Fig.  36).  These  oscillations 
are  typical  in  approximating  a  singular  field  by  a  smooth  function,  but  they  lead  to 
moderate  domain  dependence  for  the  J  integral. 

Enriching  the  EFG  trial  function  greatly  aids  in  capturing  the  singular  stress 
field  around  the  crack  tip;  when  full  enrichment  is  used  the  oscillations  are  almost 
completely  eliminated.  The  stress  along  6  =  0°  from  the  crack  tip  computed  with 
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(c)  Integration  cells 

Figure  34:  Mesh  used  for  near-tip  crack  problem  without  enrichment  using 
rings  at  r  =  0.03a,  0.09a,  0.18a. 
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Distance  from  crack  tip 


Figure  35:  Stresses  ahead  of  the  crack  tip  (6  —  0,  r  >  0)  for  the  near-tip 
crack  problem. 


Figure  36:  Angular  variation  in  stress  at  a  distance  r  =  0.1c  from  the  crack 
tip.  Solution  using  linear  basis  with  mesh  shown  in  Fig.  34. 
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Distance  from  crack  tip 


Figure  37:  Stresses  ahead  of  the  crack  tip  for  the  near-tip  crack  problem  with 
extrinsic  MLS  enrichment. 


extrinsic  MLS  enrichment  is  shown  in  Fig.  37  and  the  angular  variation  of  stress  at  r  = 
0.1a  from  the  crack  tip,  where  a  is  the  length  of  the  crack,  is  shown  in  Fig.  38.  These 
results  .were  calculated  using  the  background  nodal  arrangement  shown  in  Fig.  34a 
without  refinement  at  the  crack  tip.  The  diffraction  method  with  A  =  1  (unless 
otherwise  specified)  is  used  with  the  enriched  approximations  to  provide  smooth  and 
continuous  shape  functions  near  the  crack  tip.  It  can  be  seen  that  the  enriched  EFG 
method  is  able  to  capture  the  singularity  and  eliminate  oscillations  at  the  crack  tip 
without  using  extra  refinement  in  the  crack  tip  region.  The  stress  profiles  for  extrinsic 
and  fully  intrinsic  enrichment  are  identical  for  this  problem.  Radial  y/r  enrichment  of 
the  intrinsic  basis  captures  the  singular  stress  much  better  than  the  linear  basis  with 
no  extra  nodes,  but  it  tends  to  underestimate  the  stress  field,  as  shown  in  Fig.  39. 

The  relative  error  in  strain  energy  is  shown  in  Fig.  40,  and  it  can  be  seen  that 
the  enrichment  greatly  increases  the  absolute  accuracy.  For  a  linear  basis,  the  slope 
of  the  line  represents  the  rate  of  convergence  of  the  approximation  for  uniform  nodal 
refinement.  For  problems  with  a  singularity,  the  rate  of  convergence  is  controlled  by 
the  order  of  the  singularity  for  a  polynomial  basis  of  any  order  (Krysl  and  Belytschko, 
1996;  Strang  and  Fix,  1973).  The  rate  of  convergence  for  the  linear  basis  in  Fig.  40 
is  0.58,  which  is  just  slightly  higher  than  the  theoretical  value  of  1/2.  For  the  fully 
enriched  approximations  (extrinsic  or  intrinsic) ,  the  exact  solution  is  contained  in  the 
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Angle,  degrees 


Figure  38:  Angular  variation  of  stress  along  constant  radius  (r  =  0.1a)  from 
the  crack  tip.  Calculated  using  extrinsic  MLS  enrichment. 
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Distance  from  crack  tip 

Figure  39:  Stresses  ahead  of  the  crack  tip  (6  =  0,r  >  0)  for  the  near-tip 
crack  problem. 
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Figure  40:  Convergence  for  near-tip  crack  problem.  Diffraction  method  (A  = 
2)  used  for  enriched  basis. 


trial  function  and  according  to  the  theory  of  minimum  potential  energy,  there  should 
be  no  error  due  to  approximation.  In  this  case,  the  errors  which  arise  are  due  to 
other  effects,  such  as  quadrature  error  and  discretization  of  the  essential  boundary 
conditibns. 

Stress  intensity  factors  were  computed  using  the  domain  form  of  the  J  integral 
(Moran  and  Shih,  1987)  with  the  domain  shown  in  Fig.  41.  The  J  integral  is  theoret¬ 
ically  domain  independent;  however,  oscillations  such  as  those  apparent  in  Figs.  35 
and  36  and  inaccuracy  in  the  approximation  of  the  singular  crack-tip  stress  field  as  in 
Fig.  39  lead  to  some  domain  dependence.  Enriching  the  EFG  method  yields  domain 
independence  by  capturing  the  singularity  and  eliminating  oscillations.  Figs.  42  show 
stress  intensity  factors  plotted  against  domain  size  for  the  J  integral  with  enrichment 
by  full  and'  y/r  enrichment  versus  EFG  with  a  linear  basis  using  nodal  meshes  of  9  x  9, 
17  x  17  and  25  x  25  nodes.  Full  enrichment  shows  virtually  no  dependence  on  either 
domain  size  or  number  of  nodes  for  this  problem  and  for  y/r  enrichment,  the  stress 
intensity  factor  varies  less  than  3%  for  the  coarsest  mesh  tested. 
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Figure  41: 


Domain  used  for  evaluating  domain  form  of  J  integral. 


5.6.2  Shear  edge  crack 

In  this  example,  we  consider  a  plate  clamped  on  the  bottom  and  subjected  to  shear 
traction  r  =  1  psi  on  the  top  and  containing  an  edge  crack  of  length  a  =  Wf  2  =  3.5  in 
(see  Fig.  43).  The  material  constants  used  are  E  =  30  x  106psi  and  u  =  0.25  and 
plane  strain  conditions  are  assumed.  The  reference  solution  for  the  stress  intensity 
factors  is  (Wilson,  1969) 

k\  =  34.0  psiVm 

k2  =  4.55  psiVm . 

Stress  intensity  factors  are  computed  by  a  mixed  mode  J  integral  formulation  which 
uses  the  near  crack  tip  fields  as  auxiliary  fields;  contour  integrals  in  this  method  are 
converted  to  domain  integrals  using  Gauss’  theorem  (Moran  and  Shih,  1987). 

The  EFG  method  has  been  shown  to  work  well  for  this  problem  (Belytschko,'  Gu, 
and  Lu,  1994),  but  a  great  deal  of  nodal  refinement  near  the  crack  tip  was  required;  in 
addition  a  large  domain  was  needed  for  the  J  integral.  The  enriched  EFG  formulation 
for  this  problem  gives  results  which  are  accurate  and  require  much  less  computational 
expense.  Table  1  gives  results  for  a  linear  basis  and  enriched  EFG  using  nodal  grids 
of  57  (5  x  11)  nodes  and  324  (11  x  29)  nodes,  all  evenly  spaced.  A  domain  of  ±1  in. 
from  the  crack  tip  is  used  for  computing  the  J  integral. 
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Figure  42:  Domain  dependence  of  normalized  stress  intensity  factor  for  linear 
basis  and  enriched  EFG  for  the  near-tip  crack  problem.  Results 
are  identical  for  full  intrinsic  and  extrinsic  enrichment. 
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;  Figure  43:  Geometry  and  nodal  configuration  for  shear  edge  crack. 


Linear  Extrinsic  Intrinsic  Intrinsic 

(no  rings)  MLS  (Full)  (Vr) 

nodes  ki  k2  k\  k2  ki  k2  k\  k2 

5  x  11  26.90  3.43  34.43  4.58  32.24  4.21  31.86  4.20 

11  x  29  32.51  4.33  34.21  4.55  33.87  4.54  33.01  4.59 

Table  1:  Stress  intensity  factors  for  shear  edge  crack  calculated  by  linear 
basis  and  enriched  EFG. 
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ty=  10000  psi 


Figure  44:  Schematic  drawing  of  a  tension-loaded  plate  with  two  holes  and 
two  fatigue  cracks  growing. 


6  Numerical  Results 

6.1  Plate  with  2  Holes  and  2  Cracks 

Fatigue  crack  growth  from  rivet  holes  is  a  common  problem  for  structures  made  of 
thin,  riveted  sheets.  For  example,  an  aircraft  fuselage  is  generally  constructed  of 
aluminum  sheets  which  are  riveted  together  at  lap  joints  using  thousands  of  rivets. 
The  holes  act  as  stress  raisers  to  create  a  local  stress  concentration  and  small  defects 
often  develop  in  to  cracks.  Fig.  44  presents  a  simplified  model  of  a  plate  with  two 
holes  subjected  to  a  remote  tension  field.  Initial  cracks  are  assumed  to  emanate  from 
each  hole,  with  the  left  crack  at  an  angle  a  and  the  right  crack  at  an  angle  j3  as  shown. 
Placing  the  cracks  at  an  angle  allows  for  the  investigation  of  crack  passing  and  crack 
bridging  where  the  crack  tips  will  actually  avoid  each  other  (Melin,  1983). 

The  first  case  shown  is  for  a  =  j3  =  45°.  Figs.  45  shows  Mises  stress  contours  as 
the  cracks  are  passing  each  other  and  after  the  cracks  have  passed.  Before  the  cracks 
pass,  there  is  significant  load-bearing  capacity  in  the  region  between  the  holes  and 
up  to  the  point  where  the  cracks  pass  the  crack  tips  are  the  area  of  highest  stress. 
However,  after  the  cracks  pass,  the  region  between  the  holes  is  severely  weakened  and 
the  high  stress  area  is  shifted  to  the  edge  of  the  holes  away  from  the  cracks.  This 
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(a)  step  8  (b)  Step  16 

Figure  45:  Mises  stress  contours  for  a  plate  with  two  holes  and  two  cracks 
growing.  Initial  cracks  are  at  angles  a  =  [3  =  45°. 


is  further  illustrated  by  the  stress  intensity  factors  in  Fig.  46.  For  the  first  several 
steps  of  crack  growth,  the  mode  1  stress  intensity  factors  are  increasing,  but  after  the 
cracks  pass  the  mode  1  stress  intensity  factors  decrease.  The  mode  2  stress  intensity 
factors  are  a  measure  of  crack  smoothness  and  curvature  and  it  can  be  noted  that  as 
the  cracks  pass,  they  begin  to  curve  and  the  mode  2  stress  intensity  factor  deviates 
slightly  from  zero. 

The  next  case  shown  is  for  a  =  15 °,/3  = ‘5°.  At  this  shallower  angle  the  crack 
interaction  is  much  more  significant.  Figs.  47  shows  Mises  stress  contours  at  two 
stages  of  crack  growth.  These  results  show  a  similar  trend  as  in  the  previous  example 
where  the  crack  propagation  weakens  the  area  between  the  holes  to  the  point  where 
it  does  not  carry  load  any  longer.  It  is  interesting  to  note  that  even  at  this  shallow 
angle  the  crack  tips  do  not  coalesce,  but  rather  avoid  each  other  only  to  curve  around 
and  meet  each  other  behind  the  tip.  This  phenomena  is  known  as  crack  bridging 
(Melin,  1983). 
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Figure  46:  Stress  intensity  factors  for  plate  with  two  holes  and  two  cracks. 
Initial  cracks  are  at  angles  a  =  45°,  j3  =  —45°. 


(a)  step  7 


(b)  step  11 


Figure  47:  Mises  stress  contours  for  crack  growth  in  a  plate  with  two  holes. 
Initial  crack  angles  are  a  =  15 °,(i  =  5°. 
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Figure  48:  Stress  intensity  factors  for  two  cracks  growing  in  a  plate  with  two 
holes.  Initial  crack  angles  are  a  =  15°,  /3  =  5°. 


6.2  Double  Cantilever  Beam 

In  this  example,  a  quasi-static  crack  is  allowed  to  grow  in  a  double  cantilever  beam 
(DCB)  specimen.  The  geometry  is  shown  in  Fig.  49,  with  L  —  300  mm,  h  =  100  mm, 
a  =  138  mm,  and  P  —  100  N.  Plane  stress  linear  elastic  conditions  are  assumed  with 
Young’s  modulus  E  —  200  GPa  and  Poisson’s  ratio  u  =  0.3. 

In  addition,  a  small  perturbation  is  introduced  at  the  crack  tip  consisting  of  a  small 
kink  at  angle  dO  relative  to  the  direction  of  the  long  crack.  The  double  cantilever 
beam  exhibits  crack  path  instability,  meaning  that  the  crack  will  curve  away  from  the 
original  straight  path  after  a  small  perturbation.  Physically,  a  perturbation  may  arise 
from  an  inclusion  or  other  imperfection  in  the  material  or  original  crack  geometry. 
Sumi,  Nemat-Nasser,  and  Keer  (1985)  presented  experimental  results  demonstrating 
the  unstable  nature  of  the  crack  path.  Cotterell  (1966)  and  Cotterell  and  Rice  (1980) 
showed  through  perturbation  techniques  that  the  second  term  in  the  near-tip  asymp¬ 
totic  expansion  of  stress,  which  is  the  non-singular  stress  parallel  to  the  crack,  can  be 
used  to  identify  the  instability.  Sumi,  Nemat-Nasser,  and  Keer  (1985)  expanded  this 
analysis  to  show  situations  of  intermediate  stability  where  either  the  path  is  initially 
stable  and  becomes  unstable  or  else  the  path  is  initially  unstable  but  becomes  stable 
after  some  crack  growth.  Sumi  (1985)  performed  numerical  studies  using  finite  ele- 
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Figure  49:  Initial  geometry  of  double  cantilever  beam  (DCB)  specimen 


ments  and  found  that  by  using  a  small  perturbation  at  the  crack  tip  it  was  possible 
to  simulate  the  experimental  results;  however,  the  crack  path  was  incorporated  in  the 
initial  mesh  or  treated  by  remeshing. 

EFG  results  for  three  stages  of  crack  growth  are  shown  in  Fig.  50  for  a  linear 
basis  and  in  Fig.  50  for  a  linear  +\/r  basis.  A  perturbation  of  length  Ax  =  12  mm 
at  an  angle  d8  =  4.8°  is  placed  at  the  crack  tip  resulting  in  an  unstable  crack  path 
which  deviates  from  the  original  straight  path  and  curves  toward  the  boundary.  A 
star-shaped  array  of  nodes  is  placed  at  the  crack  tip  for  enhanced  spatial  resolution 
and  the  crack  growth  step  is  A  a  =  5  mm.  Figs.  51  and  52  show  the  effect  of  using 
a  crack  step  size  which  is  too  large.  The  predicted  path  will  contain  significant  error 
when  the  step  size  is  not  small  enough  to  resolve  the  sharp  change  in  the  path.  The 
stress  intensity  factors  for  the  linear  basis  and  linear  +  \fr  basis  are  shown  in  Fig.  53 
using  a  crack  increment  of  Aa  =  5  mm.  The  stress  intensity  factors  for  the  two  bases 
are  quite  similar,  accounting  for  the  similar  crack  path  results  shown  in  Figs.  50. 


6.3  Crack  Growth  from  a  Fillet 

This  example  shows  the  growth  of  a  crack  from  a  fillet  in  a  structural  member.  The 
set-up  to  be  modeled  is  shown  in  Fig.  54,  with  the  actual  domain  modeled  outlined 
with  a  dashed  line.  Sumi,  Yang,  and  Wang  (1995)  performed  experiments  on  a 
similar  structure  but  also  included  the  effects  of  welding  residual  stresses  and  varied 
the  bending  stiffness  of  the  structure  by  varying  the  size  of  the  bottom  I-beam.  The 
results  presented  here  are  for  a  simplified  model  in  which  the  residual  stresses  due  to 
welding  are  neglected.  In  addition,  only  the  limiting  cases  for  the  bottom  I-beam  of 


specimen  using  a  linear  basis  and  a  linear  +\/r  basis.  Crack 
growth  increment  is  A  a  =  5  mm. 
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Figure  52:  Crack  path  in  double  cantilever  beam  specimen  for  2  growth  step 
sizes  using  a  linear  +\/r  basis. 
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Figure  53:  Stress  intensity  factors  versus  crack  growth  for  linear  basis  and 
linear  +y/r  basis. 


rigid  (very  thick  beam)  and  flexible  (very  thin  beam)  are  considered.  The  material 
is  assumed  to  be  plane  strain  linear  elastic  with  Young’s  modulus  and  Poisson’s  ratio 
of  E  —  200  GPa  and  v  =  0.3,  respectively.  The  applied  load  is  P  =  1.0  N  and 
the  fillet  radius  is  p  =  20  mm.  The  initial  crack  length  is  a0  =  5  mm.  The  EFG 
formulation  is  enriched  using  the  enriched  basis  in  the  region  surrounding  the  crack 
tip  (n  =  0,  r2  =  12  mm)  and  the  enriched  fields  are  mapped  using  the  procedure  of 
Section  5.5. 

The  crack  path  evolution  shown  in  Fig.  55  is  for  the  case  in  which  the  supporting 
I-beam  is  very  thin.  In  this  case,  the  crack  curves  sharply  downward  and  propagates 
toward  the  bottom  of  the  structure.  In  contrast,  when  the  structure  is  supported  by 
a  rigid  I-beam,  the  crack  will  propagate  almost  directly  toward  the  opposite  fillet  as 
shown  in  Fig.  56. 

6.4  Beam  Under  3-Point  Bending 

The  next  example  is  crack  growth  in  a  beam  under  3-point  bending.  Three  holes  are 
located  in  a  vertical  line  offset  from  the  centerline  and  a  crack  is  seeded  at  the  bottom 
of  the  beam  and  allowed  to  grow  (see  Fig.  57).  Bittencourt,  Sousa,  Wawrzynek,  and 
Ingraffea  (1995)  presented  experimental  and  numerical  results  for  this  problem  which 
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Figure  54:  Experimental  set-up  from  Sumi,  et  al.  (1995). 


Configuration  Crack  length,  a  Crack  offset,  b 
Setup  A  1.5-  5.0 

Setup  B _ TO _  6.0 _ 

Table  2:  Crack  configurations  run  for  beam  under  3-point  bending. 


showed  that  based  on  the  location  and  length  of  the  original  crack,  the  path  of  crack 
growth  would  either  intersect  one  of  the  holes  or  pass  between  them. 

This  problem  is  solved  by  the  EFG  method  using  a  linear  basis  except  for  the  crack 
tip  region  where  the  basis  is  locally  enriched.  The  coupling  procedure  in  Section  5.4 
is  used'  with  T\  =  0.0001,  =  0.5.  The  locally  enriched  basis  has  the  advantage 

of  yielding  accurate  results  for  stress  intensity  factors  without  additional  refinement 
at  the  crack  tip.  The  presence  of  the  holes  makes  this  a  useful  feature  because  the 
additional  nodes  would  create  difficulties  if  they  fell  inside  the  holes. 

Table  2  gives  the  dimensions  and  location  of  the  initial  crack  for  the  two  cases  run 
and  Figs.  58  shows  EFG  results  corresponding  to  these  two  cases.  For  the  first  case 
solved,  the  crack  passed  between  the  bottom  two  holes  and  then  passed  very  close 
to  the  middle  hole  on  the  opposite  side.  Experimental  results  showed  that  the  crack 
actually  curved  into  the  hole;  however,  this  effect  was  not  captured  numerically  by 
either  EFG  or  by  the  finite  elements  with  remeshing  as  shown  in  Fig.  59.  For  the 
second  setup,  the  crack  was  moved  farther  from  the  centerline  and  shortened.  Results 
showed  that  the  crack  grew  directly  toward  the  center  hole  for  this  case. 

6.5  Plate  with  31  Holes 

The  examples  presented  to  this  point  have  been  for  relatively  simple  geometries. 
To  demonstrate  the  power  and  flexibility  of  EFG,  a  problem  with  a  more  compli- 
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Figure  55:  Crack  path  progression  from  a  fillet  for  the  case  of  a  thin  I-beam. 
Crack  tip  field  enriched  by  enriched  basis  for  r  =  0  —  12  mm. 
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Figure  56:  Crack  path  from  a  fillet  for  the  case  of  a  rigid  I-beam.  Crack  tip 
field  enriched  by  enriched  basis  for  r  =  0  —  12  mm. 


Figure  57:  Schematic  drawing  of  beam  with  three  holes  subjected  to  3-point 
bending. 
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(a)  Setup  A 


(b)  Setup  B 


Figure  58:  Final  crack  growth  results  for  two  different  initial  crack  configu¬ 
rations. 
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Figure  59:  Experimental  and  numerical  results  from  Bittencourt,  Sousa, 
Wawrzynek,  and  Ingraffea  (1995)  for  setup  A. 
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Figure  60:  Two  experimental  specimens  for  crack  propagation  in  a  plate  with 
31  randomly  spaced  holes  from  Al-Ostaz  and  Jasiuk  (1995). 


cated  geometry  is  solved.  Consider  a  specimen  with  thirty-one  holes  randomly  spaced 
throughout  the  domain.  Crack  growth  in  such  a  specimen  is  quite  challenging  and 
several  of  the  enhancements  for  meshless  methods  are  useful  for  solving  this  problem. 

Al-Ostaz  and  Jasiuk  (1995)  presented  experimental  results  for  brittle  fracture 
propagation  in  specimens  with  this  geometry.  The  cracks  began  at  the  edge  of  holes 
and  tended  to  intersect  the  holes  directly  until  the  whole  specimen  had  fractured 
(see  Figs.  60).  They  also  performed  finite  element  studies  with  a  fine  mesh  in  which 
fracture  was  simulated  by  the  removal  of  elements.  Whenever  the  strain  energy  or 
the  maximum  principal  stress  in  an  element  reached  a  critical  value,  the  element  was 
deleted  from  the  mesh,  thereby  extending  the  crack. 

This  problem  was  modeled  by  EFG  using  an  initial  nodal  distribution  of  2396 
nodes.  A  stress  analysis  was  performed  with  no  cracks  in  order  to  determine  where  the 
highest  stresses  were.  Small  cracks  were  seeded  at  the  seven  highest  stress  locations. 
The  bottom  edge  was  fixed  and  a  constant  displacement  was  applied  in  the  vertical 
direction  to  the  top  surface.  The  displacement  along  the  top  surface  was  incremented 
each  step  by  0.02%  each  step.  The  enriched  basis  from  Eq.  (50)  was  used  in  the  region 
surrounding  each  crack  and  was  ramped  to  a  linear  basis.  The  enriched  basis  was 
chosen  over  increased  spatial  refinement  because  of  the  complex  geometry. 

Results  following  12  steps  of  crack  growth  are  shown  in  Fig.  62.  After  a  few 
steps  of  growth,  a  few  dominant  cracks  developed  and  the  other  cracks  arrested. 
Some  cracks  grew  into  other  holes,  but  in  general  the  cracks  seemed  to  deflect  away 
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Figure  61:  Mises  stress  contours  for  stress  analysis  without  cracks  for  plate 
with  31  randomly  spaced  holes.  Cracks  are  seeded  at  the  seven 
highest  stress  locations. 


from  the  holes.  This  result  is  in  contrast  to  the  experimental  and  numerical  results 
of  Al-Ostaz  and  Jasiuk  (1995)  who  found  that  the  crack  tended  to  grow  directly 
into  the  holes.  Possible  explanations  of  the  differences  with  the  experiments  and 
the  EFG  computations  is  dynamic  effects.  The  crack  speed  in  the  experiments  was 
quite  high  (the  specimens  fractured  completely  in  less  than  1/3000  sec)  while  the 
EFG  computations  assumed  quasi-static,  crack  propagation.  The  cracking  patterns 
obtained  by  EFG  use  the  principles  of  fracture  mechanics,  whereas  the  finite  element 
cracking  pattern  was  predicted  by  the  area  of  high  stress. 


7  Conclusions 

A  meshless  method  called  the  element-free  Galerkin  (EFG)  method  has  been  pre¬ 
sented  as  a  computational  tool  for  simulating  fatigue  and  quasi-static  crack  propaga¬ 
tion.  EFG  approximations  are  formed  from  a  basis,  polynomial  or  enriched,  and  a  set 
of  unknown  coefficients.  A  moving  least  squares  (MLS)  method  is  used  to  determine 
the  unknown  coefficients  by  minimizing  a  weighted  £2  norm  written  in  terms  of  nodal 
coefficients.  The  £2  norm  in  multiplied  by  a  weight  function  with  compact  support, 
resulting  in  a  local  approximation,  i.e.,  the  approximation  is  written  in  terms  of  the 
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(a)  Initial  setup 


(b)  After  12  steps 


(c)  Mises  stress  contours 


Figure  62: 


Quasi-static  crack  propagation  in  plate  with  31  randomly  spaced 
holes. 
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surrounding  nodes  so  the  discrete  equations  are  sparse.  A  weight  function  is  defined 
at  each  node  and  common  support  shapes  are  circles  and  squares.  A  crack  is  modeled 
as  a  line  segment  in  the  interior  of  the  domain.  When  a  domain  of  influence  of  a  node 
near  the  crack  intersects  the  crack,  it  is  truncated  along  the  crack  surface.  Crack 
growth  is  modeled  by  adding  an  additional  segment  at  the  tip. 

Smoothing  methods  for  meshless  approximations  near  nonconvex  boundaries  are 
presented.  Boundaries  in  meshless  methods  are  enforced  by  truncating  the  domains 
of  influence  of  the  nodal  weight  functions  along  the  boundary.  For  nonconvex  bound¬ 
aries,  this  leads  to  discontinuities  which  extend  into  the  domain.  For  the  case  of  a 
crack,  the  weight  function  of  a  node  near  the  crack  tip  will  be  discontinuous  along  the 
ray  from  the  node  which  grazes  the  crack  tip.  Because  the  shape  functions  consist 
of  the  contributions  of  several  weight  functions,  all  shape  functions  near  the  crack 
tip  will  contain  discontinuities  along  this  line.  These  discontinuities  are  dealt  with  in 
several  ways.  The  first  is  to  simply  ignore  them  and  use  the  discontinuous  approxima¬ 
tions.  The  second  technique  is  to  use  the  diffraction  method  in  which  nodal  domains 
of  influence  for  nodes  near  a  crack  tip  or  other  nonconvex  boundary  are  wrapped 
around  the  boundary  similar  to  the  way  in  which  light  diffracts.  A  third  method, 
called  the  transparency  method,  is  developed  only  for  cracks.  In  this  method,  the 
domains  of  influence  for  nodes  near  the  crack  tip  are  not  severed  at  the  crack  tip,  but 
are  truncated  gradually  until  at  a  short  distance  behind  the  tip  they  are  truncated 
completely.  Another  method  for  obtaining  continuous  approximations  is  called  the 
“see-through”  method  or  the  continuous  line  criterion.  According  to  this  criterion,  if 
a  continuous  line  can  be  drawn  between  two  points  without  leaving  the  domain  of  in¬ 
fluence,,  then  the  weight  function  is  not  modified.  This  method  works  well  for  smooth 
boundaries,  such  as  holes,  for  refined  nodal  discretizations,  but  is  not  recommended 
for  sharp  boundaries  such  as  cracks. 

Techniques  for  enriching  meshless  approximations  for  linear  elastic  fracture  are 
presented.  Enrichment  can  be  used  when  the  form  of  the  solution  is  known;  this  in¬ 
formation  may  be  used  to  form  better  approximations  yielding  better  solutions  with 
fewer  degrees  of  freedom.  The  enrichment  is  carried  out  in  two  basic  ways,  extrinsi- 
cally  and  intrinsically.  Extrinsic  enrichment  involves  adding  the  form  of  the  solution 
to  the  EFG  approximations  constructed  from  a  standard  polynomial  basis.  Within 
extrinsic  enrichment,  extrinsic  MLS  enrichment  and  extrinsic  PU  enrichment  are  pre¬ 
sented.  Extrinsic  MLS  enrichment  consists  of  adding  the  enrichment  functions  to 
the  approximation  which  is  used  in  the  MLS  methodology.  Extrinsic  PU  enrichment 
consists  of  using  the  EFG  shape  functions  as  partitions  of  unity  (hence  the  term  PU) 
and  adding  the  enrichment  terms  in  an  extrinsic  basis  which  is  multiplied  by  the 
partition  of  unity  and  some  unknown  nodal  coefficients.  The  other  form  of  enrich¬ 
ment,  which  is  called  intrinsic  basis  enrichment,  consists  of  expanding  the  existing 
polynomial  basis  to  include  the  terms  from  the  asymptotic  near  tip  solution.  Both 
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full  and  partial  enrichment  are  presented.  Full  enrichment  adds  the  complete  span  of 
the  near  tip  displacement  field  to  the  intrinsic  basis,  while  partial  enrichment  adds 
only  the  y/r  radial  dependence  to  the  basis.  Computing  enriched  approximations  is 
more  expensive  and  is  unnecessary  away  from  a  crack  tip.  To  this  end,  a  consistent 
coupling  procedure  is  defined  which  allows  enriched  approximations  to  be  used  near 
a  crack  tip  and  then  blended  with  the  approximations  from  a  linear  basis.  The  cou¬ 
pling  is  defined  such  that  consistency  is  maintained  and  no  discontinuities  arise  in 
the  displacement  field.  When  the  crack  is  kinked,  a  transformation  is  defined  which 
will  map  the  enriched  fields  to  coincide  with  the  actual  crack  path.  This  modification 
is  required  because  the  enriched  field  includes  a  discontinuity  at  6  =  ±7r  which  does 
not  coincide  with  the  required  discontinuity  beyond  the  kink. 

The  element-free  Galerkin  method  is  shown  to  be  an  effective  method  for  solv¬ 
ing  solid  mechanics  problems.  The  method  is  easily  able  to  handle  arbitrary  crack 
propagation. 
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